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EXECUTIVE SUMMARY 


This report presents the results of a feasibility study of the Quadrilateralized 
Spherical Cube Earth Data Base. The structure of this data base is derived 
from a spherical cube, which is obtained by radially projecting a cube onto its 
circumscribing sphere. An appropriate set of curvilinear coordinates is chosen 
such that the resolution cells on the spherical cube are of equal area and are 
also of essentially the same shape. This geometrical construction eliminates 
both the necessity to account for nonequal-area cells, and the need to compute 
the location of every cell in the data base. 

The main properties of the Earth data base are as follows: 

1 . The indexing scheme is binary in nature, and telescopic in the sense 
that each additional level of resolution is addressed by appending 
additional binary bits. Thus, minimal work is needed for indexing 
cells of higher resolution. 

2. The resolution cells are strung together in a two-dimensional man- 
ner, so as to accomplish area coverage with a serial bit string. 
Consequently, a higher degree of proximity is achieved for near- 
neighbors in this stringing pattern than in the usual one-dimensional 
array of stringing by rows and columns. 

3. The cell addresses are readily computed because of the indexing 
scheme which is the same regardless of the resolution level, and 
because of the stringing pattern which permits the decomposition 
of the cell address into two independent binary indices. 

4. The conversion from geographic coordinates to data base coordin- 
ates is comparatively simple because of the simplicity of the data 
base structure. 
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5. Incoming data can be stored rapidly by interpolation, using bench- 
marks only occasionally. This method of fast-filling is made pos- 
sible by the equal-area nature and translational shape invariance 
of the data base resolution cells. 

6. Input/output operations with this data base are also simplified 
because of the rectangularized nature of the data base records and 
the rhombic nature of the interpolation blocks. ‘ 

7. The user can rapidly and directly access data corresponding to 
specified geographic regions of arbitrary shape and size. This 
data-accessing is accomplished by retrieving the relatively few 
bit-strings which lie within the associated data base records. The 
rapidity and directness of data access are the result of equal-area 
resolution, translational invariance, indexing scheme, stringing 
pattern, and relatively simple coordinate transformation. 

8. The primary contemplated uses of the retrieved data are mathe- 
matical computations and visual display. For the former, the 
equal-area resolution property eliminates the need to distinguish 
between density measurements and integrated measurements. For 
the latter, the quadrilateralized nature of the resolution cells on 
the spherical cube and the comparative simplicity of coordinate 
transformation both simplify and minimize the internal operations. 


Based on numerical results obtained in this study, it is concluded that imple- 
mentation of the proposed Quadrilateralized Spherical Cube Earth Data Base 
is feasible on current large-scale computer systems. 
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SECTION 1 - INTRODUCTION 


In the numerous satellites presently orbiting the Earth, enormous amounts of 
data are continuously taken of the Earth's surface and atmosphere. These 
data are of a varied nature: topography, crop distribution, sea surface tem- 
perature, cloud coverage, etc. The measurements are used by research and 
applications personnel of diverse scientific disciplines. These users usually 
employ and Earth-oriented coordinate system, such as the traditional geo- 
graphic frame of reference. Thus, it is not surprising that almost all existing 
Earth data bases have been constructed with latitudes and longitudes as grid- 
lines, either in a patched-up partial fashion or in the entire outlay. 

However, what is convenient to the user is not necessarily also efficient from 
the standpoint of data management and data processing by the computer. Effi- 
ciency is especially important because of the large amounts of data rapidly 
acquired in global coverage, the necessity to update data continually for opera- 
tional use, and the desire to access directly relatively small amounts of data 
corresponding to selected geographic regions at appropriate times. 

The high computer overhead encountered in processing can therefore be mini- 
mized by designing an Earth data base structure with constant (but selectable) 
geometric resolution cells, which are also locally invariant in shape along a 
translation in any direction. This would eliminate the necessity to account 
for nonequal-area resolution cells, and also the need to compute the location 
of every resolution cell in the data base. Moreover, the design should also 
utilize a fairly simple transformation between the user-preferred geographic 
coordinates and the internal data base coordinates. This would greatly facili- 
tate arithmetic and transfer operations desired by the user in mathematical 
computations or in graphic display. 
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In the previous study (Reference 1), six data base structures were conceptu- 
alized. Of these, the most promising was the Quadrilateralized Spherical 
Cube Data Base. In this model, the sphere is visualized as a spherical cube, 
as illustrated in Figure 1-1. This spherical cube is obtained by radially pro- 
jecting the edges of an inscribed cube, as shown in Figure 1-2. 

From Figure 1-2, it is obvious that equal-area elements on the plane square 
do not radially project as equal-area elements on the spherical square. For 
example, those elements near the center of the plane square have larger pro- 
jections than those elements near the edges of the plane square. Hence, if a 
rectangular grid of equal-area elements is first constructed on the plane 
square, it is then necessary to distort this grid into a curvilinear network so 
that the elements near the center are smaller than those near the edges. The 
distortion is such that when the curvilinear elements are projected radially, 
equal-area elements are again obtained on the spherical square. The desired 
sequence of transformations is illustrated in Figure 1-3 through 1-5. The 
mathematical details of deriving these transformations are discussed in 
Section 2. 

For the present, it suffices to say that it is possible to obtain a world map 
such as Figure 1-6. This map illustrates the continental outlines as they 
would appear on the cube with the original undistorted rectangular coordinates. 
This is accomplished by reversing the sequence of transformations previously 
illustrated by Figures 1-3 through 1-5. Thus, in Figure 1-6, equal-area 
regions correspond to equal-area regions on the spherical Earth. An examin- 
ation of this planar equal-area world map shows that the distortion of the con- 
tinental outlines is not as great as might be expected. 
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Figure 1-3. Plane Square With Cartesian Coordinates 
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Figure 1-4. Plane Square With Curvilinear Coordinates 



Figure 1-5. Spherical Square With Curvilinear Coordinates 
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SECTION 2 - MATHEMATICAL FORMULATION OF 
DATA BASE STRUCTURE 


2.1 SUMMARY 

This section is concerned with the derivation of a mapping function which leads 
to the division of a spherical cube into quadrilaterals of equal area. Specifi- 
cally, the method involves: (1) dividing a plane square into rectangles of equal 
area; (2) transforming the rectangular coordinate system on the plane square 
to a curvilinear coordinate system on the same plane square; and (3) projecting 
the new nonequal-area elements on the original plane square to a subtending 
spherical surface such that equal-area spherical surface elements are again 
obtained. In addition to the direct transformation function described, analysis 
is also performed for obtaining the inverse transformation function which re- 
verses the above process. ' 

2. 2 DERIVATION OF DIRECT MAPPING FUNCTION 

First, consider a plane surface subtended by a spherical surface with 

radius R . Let r be the vector from the center of the sphere to the given 
o 

plane. As shown in Figure 2-1, let dA be an area element on this plane, 

P 

and let 'r‘ be the vector from the center of the sphere to the area element 

dA . 

P 

Let dA be the spherical area element obtained by projecting dA radially 
s P 

onto the sphere. Then, it can be readily shown that the following relation 

between dA and dA holds: 

P s 


dA 

s 


2 3 . 

R cos (r, r ) 
o 


2 

r 

o 



( 2 - 1 ) 


where fr% ) denotes the angle between r and 

O 
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Let (£ , r? , r ) denote the components of the vector r . Then, it follows 
^ o 

that 


cos (r, r ) 
o 


/ 2 .2 2 

r o+ | 


172 " 


(2-2) 


Moreover, for convenience, let the unit of length be chosen such that the 
radius, R , of the sphere is equal to unity. Then, Equations (2-1) and (2-2) 
yield 


dA 

s 


r 

o 


+ t 




(2-3) 


Next, consider a cube together with a circumscribing spherical surface. On 
each of the six plane faces of the cube, a rectangular coordinate system (x,y) 
may be defined, the domain of definition being -r Q ^ x , y ^ r^ . It may be 
easily verified that 


1 

r = "7= 

0 x/iT 


(2-4) 


Let a new coordinate system (£ , 77 ) be defined by 

£ = f. (x, y) 
r? = T 7 (x,y) 


(2-5) 


where 4 (x , y) and r? (x , y) are independent arbitrary functions. 
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The new area element d£dT7 is related to the original area element dxdy by 


d^dr? = J dxdy 


( 2 - 6 ) 


where 



is the Jacobian of transformation 


J 



H H 

9x Sy 
5r? hr] 
Sx Sy 


(2-7) 


If this new area element is projected radially onto the surface of the sphere, 
Equations (2-3) and (2-6) yield 


dA 

s 


r 

o 


r + 4 




(2-8) 


which relates the spherical area element dA to the original area element 

s 

dxdy . For original equal-area plane elements dxdy to transform into equal- 
area spherical elements dA^ , it follows that 


r 

o 






(2-9) 
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or 


r 

o 




5? I _ ^] = X 2 

ay ay ax/ 


( 2 - 10 ) 


where X is a constant. 

2 

It is easy to verify that the value of X is equal to the ratio of the area of the 
spherical square to the area of the plane square, i.e. , 


X 


2 


2tt/ 3 
4 r 2 


n_ 

2 


( 2 - 11 ) 


An alternative form of Equation (2-10) is 


H sr? _ a|_ *n = 2 L + ii 

ax ay ay ax 'I ^2 



(2-12) 


where 


2 

y 


2 

= r X 
o 


•n_ 

6 


(2-13) 


2 

From Equations (2-7) and (2-12), it is seen that y may be interpreted as the 
area-scale of transformation at the point (£ = 0 , 77 = 0) . 
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Equation (2-12) in itself is quite general. It is now desirable to specify the 
following general properties for the transformation from (x , y) to (£ , r\) . 

1. To preserve symmetry in the transformation Equation (2-5), it is 
required that 


i = f(x, y) 
V =f(y, x) 


(2-14) 


Equation (2-14) states that £ and r) have exactly the same form of dependence 
on x and y , except that the roles of x and y are interchanged. Moreover, 
symmetry preservation also requires that the function f(x , y) be odd in x 
and even in y , i. e. , 


f(-x, y) = -f(x, y) 
f(x, -y) = f(x, y) 


(2-15) 


As a consequence of Equations (2-14) and (2-15), it is seen that the origin maps 
back into itself, i. e. , 


| f(0, y) - 0 (2-16) 

[ 2. To map points on the sides of the square back into points on the 

same sides, it is necessary that 

I 

I 

f( r » y) = r (2-17) 

k o o 
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As a consequence of all the above requirements, it may be shown that 

are zero at the points (0,0) and (r^ j r Q ) ■ Therefore, from Equa- 
tion (2-12), it follows that 


and 


K 

Sx 


x=0 

y=o 


_ St? 

sy 


x=0 

y=o 


ITT 


y=Jj= 0- 72360 


12545 582 


(2-18) 


Sx 


x=r 

o 




sy 


x=r 




6494 54166 187 


(2-19) 


If f(x , y) can be expanded in a power series in x and y , then Equation (2-16) 
requires that 


ou oo 

y)-*E Z 


i=0 j=0 


2i 2j 
a. . x y 
ij 


( 2 - 20 ) 


The condition in Equation (2-18) yields 


a 


00 


= 7 


(2-21) 


The condition in Equation (2-17) may be incorporated into f(x , y) by writing 
it in the form 


f(x, y) = yx + ^T^-x 3 + L? _ x 2 \ x ^ b x l y 2 ^ (2-22) 

r (i+j)sl 
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It may be shown that Equation (2-12), together with the conditions given by 
Equations (2-16) through (2-18), are not sufficient to determine uniquely the 
transformation in Equation (2-14). This nonuniqueness is manifested by the 
fact that there are more unknowns (b„) than equations when Equation (2-22) 
is substituted into Equation (2-12) and terms of the same degree are equated. 

Finally, to incorporate the condition in Equation (2-19), it is most efficient to 
express f(x , y) in the following form. The details for arriving at this form 
are given in Appendix A. 


f(x, y) = yx + x 3 


2/2 2 ' 
+ xy (r - x 


, / 2 2\ 2i 2j 

6 + r - y . c.. x y 

\ o /no i] 


j^o 


(2-23) 


+ x 3 (r 3 -x 2 ) « + (r 2 -x 2 ) £d 

iso 


2i 


x 


where 



= 0.79048 64491 208 


1 / 4 ' 

w= __ 3 - 2y- M - 2r 6 

2r ' 
o 

= - 1.2254 41487 984 


(2-24) 


(2-25) 
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An approximate mapping function may be obtained by truncating the series 

expansion in Equation (2-23) at some degree, and then obtaining the coefficients 

c.. and d. which minimize the following residual function: 
ij i 



fei 




\§X 

ay ' 

ay 

ax/ 


dxdy 


(2-26) 


This residual function is obtained by considering Equation (2-10) or (2-12). 
Then, 0(c.. , d.) is evidently equal to zero for the exact transformation func- 
tion f(x , y) . For computational purposes, Equation (2-26) is replaced by 



-,>-EE 


x k y i 



( .2 2\-3/2 

1 + S + ^j 


X 


/a| m\_x 2 

\9x 3y 9y 9x / 


2 


(2-27) 


where the points (x^ , y^) are chosen to form a regular grid over the plane 
square. A computer software program for performing this minimization prob- 
lem is given in Appendix E. For a second-degree approximation of the series 
in Equation (2-23), the following values of c„ and d. are obtained: 

c Q0 =-2.7217 05366 1814 


c^ -5.5842 16830 5430 
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= 2.1711 17480 9423 


C 2Q =-3.4578 62747 3390 
c =-6.4160 15152 6783 
c Q2 = 1.9736 26575 8872 
d = 1.4833 12929 4187 
d = 1.1199 72606 9742 
d = 6.0515 38216 1464 


The corresponding mapping function f(x , y) is accurate to about five signifi- 
cant figures. 

2.3 DERIVATION OF INVERSE MAPPING FUNCTION 

Corresponding to the symmetrical direct mapping function expressed in Equa- 
tion (2-14), it may be verified that the inverse mapping function is also sym- 
metrical, i.e., 


x = f* (£, T?) 
y = f* (n, £) 


(2-28) 
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As discussed in Appendix A, f* (£ , rj) must be expressed in the form 


f* (|, 17 ) = y*Z + 


(1 - V*) 

2 

r 

o 


(2-29) 


»■ « - <) 



jasO • 


-* + ( r ° - <*) S v ?2 ‘] 


where 



y * - 

r l 


1 

4 . 
y + r 0 
o 


* - . 


1 _ 4 . 

u+ 2r 0 
o 


5* 


(M 1 + - M*) 
2r 4 


CO 51 


= (3 - 2y* ~ fj.* - 2r 4 6*j 

2r 


6 * = — 
1 2 


'(y x * - y*) 


- 6 * 


(2-30) 
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An approximate inverse mapping function may be obtained by truncating the 
series expansion in Equation (2-29) at some degree, and then obtaining the 
coefficients c * and d * which minimize the following residual function: 

ij 1 


0 *(c..*, d*)= f°f ° Jj^x - f* (f(x,y), f(y,x)) 
r o r o 

\l/2 

+ jy - f* (f(y,x), f(x, y))j 2 j dx dy 


(2-31) 


In obtaining this residual function, the direct mapping function f(x , y) is con- 
sidered, as given by Equation (2-23). Then, 0* (c..* , d.*) is evidently equal 
to zero for the exact inverse mapping function f* (£ , 77) . Again, for compu- 
tational purposes, Equation (2-31) is replace by 


0*(C 


iJ 


d.*) = ^ £ |[x-f* (f(x, y), f(y,x)) 

\ y x ( 


+ 


y - f* (f(y. 


x), 



(2-32) 


A computer software program for performing this minimization problem is also 
given in Appendix E. For a second-degree approximation of the series in Equa- 
tion (2-29), the following values of c„* and d^* were obtained: 

c * - 3. 973 89249 

c 1Q *= 6.591 19476 

C 01* = ~ 25 ' 368 92536 

c * = -73.064 97000 

Ci U 
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c *= 77.381 61133 


° 02 * = 21,685 89623 

d * = 1.811 28250 
0 

d * = 37.635 47857 

d * = 63.000 23655 
z 

The corresponding mapping function f* (£ , r?) is accurate to about five signi- 
ficant figures. 
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SECTION 3 - ORGANIZATION OF DATA BASE 


3.1 SUMMARY 

This section mainly is concerned with the organization of the data base, whose 
basic structure is derived from the equal-area division of the spherical cube. 

It suffices, therefore, to limit most of the discussion to just one spherical 
square. The binary indexing scheme of the data base will be described. The 
reasons for adopting this indexing scheme are then discussed. This is then 
followed by an outline for computing the serial string index for the cell 
addresses. 

3.2 INDEXING OF DATA BASE 

The nature of the data base system adopted is extremely simple since all data 
are stored in a single large array. This organization is natural since: 

1. The contemplated data base corresponds to a completely filled 
array 

2. No data item has other than a positional relationship to another 

Thus, all data in the data base may be accessed by directly calculated positional 
indices with no storage inefficiencies due to unused areas of the array. 

The mapping of the data base onto storage may be accomplished in a number of 
ways. The classic method is mapping by columns and rows, as is done in 
FORTRAN array storage and as implied in the usual mathematical matrix 
notation. Other mapping techniques include scatter storage in which a psuedo- 
random, one-to-one mapping is employed as well as various consistent logical 
stringing techniques. 

The specific scheme studied in this report is based on the binary division 
scheme discussed in Reference 1. This mapping starts at the level of the faces 
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of the spherical cube, numbering these faces 1 through 6 as in Figure 3-1. 

Each face is divided, to the requisite resolution level, by a two-dimensional 
binary grid, as shown in Figure 3-2. On each level of division, the areas are 
divide^ into quadrants, which are labeled by a 2-bit binary number. Each 
level of division, n , is indicated by the addition of two binary bits to the least- 
significant end of a 2N bit binary number. This binary index is the serial 
location of a point in the 2N location data array. 

An address to the third level of division would appear as in Figure 3-3. 

This organization will be referred to as the serial addressing scheme, and the 
data array as the serial data string. 

3.3 REASONS FOR ADOPTING INDEXING SCHEME 

In developing a method of ordering data in the data base, there are three 
obvious advantages to using a serial addressing scheme rather than a normal 
column and row addressing scheme: 

1. Reduction in I/O time through maintenance of near-neighbor rela- 
tionships 

2. Compactness of arrays containing addresses 

3. Maintenance of a consistent addressing scheme regardless of reso- 
lution level 

The serial addressing scheme reduces I/O time for disk type storage devices 
because more near neighbors of a point are within the range which requires no 
arm motion for accessing. A specific example is provided in Section 3.5. 

The expression of addresses as a single bit string allows storage of addresses 
as single machine words, whereas a two-dimensional addressing scheme would 
require two or three words, including one for the face number. Finally, the 
expandibility and generality of the serial string permits the use at any resolution 
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Figure 3-1. Face Numbering Scheme 
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Figure 3-3. Illustrative Labeling by Binary Bits 
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level without regard to physical storage considerations, such as record size. 
Any reasonable matrix type storage scheme would require a dual (or multiple) 
level of addresses for record and item within record location in the serial 
scheme. This is accomplished simply by considering the high order m bits 
as the record number, and the low order n - m bits as the address within 
record. 

A major advantage discovered while studying implementation of this scheme was 
the extreme rapidity of index calculation. This last advantage is gained through 
changing the basic linking from the previously proposed pattern, shown in 
Figure 3-4, to that shown in Figure 3-5. This pattern retains all of the advan- 
tages of the former scheme, and also permits a simpler calculation of the 
serial index. 

The motivation for this pattern may be best understood by considering the data 
base as a space filling binary tree structure in two dimensions. Figure 3-6 
illustrates a simple one-dimensional tree and an associated binary numbering 
scheme. Figure 3-7 extends this system to two dimensions, from which it is 
seen that it is possible to combine the two coordinates into a one-dimensional 
index. This is accomplished by utilizing the even bits of a binary number to 
represent one coordinate and the odd bits to represent the other. Table 3-1 
illustrates the values assumed by each coordinate. 

For example, the addresses of the cells (1, 1) and (3, 5) will be represented 
in binary notation as follows: 

( 1 , 1 ) - 01 + 10 = 11 

(3, 5) - 101 + 100010 = 100111 

The x and y coordinates may be recovered from the index by simply masking 
out the even or odd bits. 


3-5 



1 4 



2 3 


Figure 3-4. Previously Proposed Stringing Pattern 
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Figure 3-5. Presently Proposed Stringing Pattern 
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Table 3-1 


Binary Representation of Coordinates 


DECIMAL VALUE 

BINARY 

X-COORDINATE 

(ODD) 

Y-COORDINATE 

(EVEN) 

1 

1 

01 

10 

2 

10 

100 

1000 

3 

11 

101 

1010 

4 

100 

10000 

100000 

5 

101 

10001 

100010 

6 

110 

10100 

101000 

7 

111 

10101 

101010 
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In summary, the binary indexing scheme and the bit-string ordering for the 
first two levels of division are shown in Figure 3-8 and 3-9. 

This linkage yields the desired serial string at any level. It also allows a com- 
putation of the serial address which requires only two table look-ups and an 
add operation. 

3.4 CALCULATION OF SERIAL STRING INDEX 

As mentioned in the previous section, the serial index is the composite of the 
x and y indices. The most rapid way of obtaining these indices is by table 
look-up. The construction of the y {even bit) table may be done as follows: 

1. Looping over all bits in the table index 

2. Examining the index bit by bit (when bit n in the index is on, set 
the bit 2n in the table entry) 

The x (odd bit) table may be obtained from the even bit table by right-shifting 
the even bit table entries (Table 3-1). 

The calculation of the serial string index may also be accomplished by the con- 
struction of a very simple hardware device. This device would consist of 
three registers: an x register, a y register, and a serial string or s 
register. 

Two register-to-register instructions would provide packing from x , y to 
s and unpacking s to x , y . These instructions would initiate parallel trans- 
fer from the two n-bit coordinate registers to the 2n-bit serial register and 
vice versa. The interconnection is shown in Figure 3-10. 

3.5 NEAR-NEIGHBOR DISTANCES IN SERIAL STRING 

The serial stringing described in Section 3.3 leads to a reduction in disk I/O 
time because of a reduction of the number of seek operations. This happens 
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because the number of large gaps separating typical adjoining records is fewer 
for the serial string than for a matrix type storage scheme. 

Figure 3-11 illustrates a resolution element and its eight adjoining spatial 
neighbors. The circle of radius, r = a , is the locus of points one resolution 
element away from the current point. For this situation, the probabilities of 
the neighbor being horizontally, vertically, or diagonally adjacent are 

P- = 1/3 
h 

P_ = 1/3 
v 

p a ■ 1/3 

The individual probabilities for a specific neighbor are 

P h = 1/5 

P = 1/6 
v 

P d - 1/6 

The matrix storage scheme (column and row storage) maintains a near-neighbor 
relationship only along the horizontal or vertical axis. Figure 3-12 is a graph 
of the probability of finding a near neighbor at a given serial separation in an 
n x n matrix arranged in columns and rows. The same graph for the serial 
string described in Section 3.3 (for a 64 x 64 array) appears in Figure 3-13. 
Depending on the definition of nearness in the serial string and the size of the 
array, a gain of 50 percent to 100 percent can be realized by using the serial 
string scheme. 

The case of multi-element record stringing differs somewhat from the simple 
resolution elements relationship described above. In propagating the 
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Figure 3-11. Resolution Element and Near Neighbors 
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Figure 3-12. Separation of Near Neighbors in Matrix Storage 
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Figure 3-13. Separation of Near Neighbors in Serial String 
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interpolation vector, r is much less than a , and P^ becomes very small. 

P- and P_ remain equal and are approximately equal to l/2. Taking ±8 rec- 
ti v 

ords as a reasonable cylinder size, an advantage of 7/5 is found for the serial 
scheme. 
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SECTION 4 - DATA STORAGE 


4.1 SUMMARY 

This section is concerned with storing data into the data base by the method of 
"fast-filling. " The method involves the determination of benchmark points 
corresponding to specific scan resolution elements (pixels). The filling of 
neighboring data base points corresponding to neighboring pixels is then achieved 
by interpolation, using appropriate combinations of data base vectors and scan 
element vectors. This method of data base filling is possible because of the 
local translational invariance of the division of the spherical cube described 
in Section 2. 

4.2 CONVERSION FROM GEOGRAPHIC COORDINATES TO DATA BASE 
COORDINATES 

The basic objective of the fast-filling method is the insertion of data into the 
data base with the minimum of computation. For this purpose, it is assumed 
that over some small region the mapping of scan lines onto the data base yields 
nearly straight lines, and that the location of points along these lines can be 
found by interpolation. It is also assumed that there is no editing of data base 
information. Rather, all data are inserted into the data base as though it were 
initially empty. 

For defining the ground location of a given pixel, the following three basic sets 
of information are necessary: 

• Satellite orbital position 

• Satellite attitude 

• Time of scan 

The determination of ground location (or geographic coordinates) as a function 
of the above parameters is described in Appendix B. Once the geographic coor- 
dinates are available, it is straightforward to obtain the corresponding set of 
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Cartesian coordinates (a , b , c) in the three-dimensional space. The subse- 
quent mapping of this point into a point onto the cube representing the data base 
requires the following four operations: 

1. Determination of the face onto which the point is to be mapped 

2. Projection from the Cartesian (a , b , c) system into the curvilinear 
(£ » V) system on a plane square 

3. Transformation from (4 , V) into Cartesian (x , y) coordinates on 
a plane square 

4. Conversion from face number and (x , y) to serial location in the 
data base 

A simple method of face determination arises naturally when a suitable rela- 
tionship is chosen between the Cartesian system and the data base coordinate 
system. This system, illustrated in Figure 4-1, is one where each face of the 
data base cube is perpendicular to a major axis of the Cartesian system whose 
origin is at the center of the cube. 

The face to which a point (a , b , c) belongs is determined by finding the com- 
ponent with the largest absolute value. This point is then projected into the 
face containing the axis corresponding to the largest component and to its sign. 
It should be noted at this point, however, that a point can be projected onto any 
of three planes, two of which are extensions of other faces. This projection 
is illustrated in Figure 4-2. These projections onto extended faces will have 
application later. 

The projection from the (a , b , c) system onto the face of the cube requires the 
adoption of a consistent set of face coordinates. The set of face coordinate and 
face numbering used in Section 3 and their relation to the Cartesian system is 
illustrated in Figure 4-3. The projection is illustrated in Figure 4-4. 
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From the geometry in Figure 4-4, it is seen that the cube centered Cartesian 
components of the point (£' , 77 f , P) project onto the cube face (4 , rj) as 
follows: 


4 =4' p 
7] ~ T}' P 


(4-1) 


where P = r /p 
o 

From Figure 4-3, Table 4-1 is constructed illustrating relations between the 
(a , b , c) and (4' , 77 ’ , p) systems. 


Table 4-1. Relations Between (a , b , c) and 
(4' , r? ' , p) Systems 


FACE NUMBER 

r 

V 

P 

1 

b 

C 

a 

2 

-b 

c 

—a 

3 

—a 

c 

b 

4 

a 

c 

-b 

5 

—a 

-b 

c 

6 

—a 

b 

— c 


Using these relations, any point (a , b , c) may be mapped into an appropriate 
(4 , V) system. 
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The transformation from the curvilinear (£ , rj) system to the rectangular data 
base coordinates (x , y) has been described in Section 2. The relevant equa- 
tions are 


x = f*(£ , V) 
y = f*(rj , 4) 


(4-2) 


The serial data base location is obtained by table look-up as described in Sec- 
tion 3 or by direct calculation. 

Clearly, the application of this mapping procedure to every data point would be 
inefficient. Consequently, it is desirable to seek methods for reducing the 
number of computations. 

4.3 METHOD OF INTERPOLATION 

The simplest method for reducing the number of computations is to utilize an 
interpolation method. This interpolation may be done at any one of a number 
of stages of the mapping from scanner coordinates to data base coordinates. 

The greatest reduction in computation is achieved by performing, in the data 
base coordinate system, a linear interpolation both along the scan line and 
across scan lines. Neglecting edge effects, interpolation in this frame reduces 
the number of transformation by a factor N x M , where N is the number of 
pixels between calculated points, and M is the number of scan lines between 
calculated points. These points will be referred to as benchmarks, and the 
grid they form as the benchmark grid. 

Interpolation (Figure 4-5) is accomplished by taking three points on the bench- 
mark grid and constructing displacement vectors in the along-scan and across- 
scan directions. 
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Figure 4-5. Interpolation Using Displacement Vectors 
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The unit displacement vectors a , b are obtained from 
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IT = "b/m 


(4-3) 


Consider a pixel at point P offset by n pixels and m scan lines from the 
benchmark O . Then its coordinates are given by 

OP = ha“+ mb (4-4) 


Two factors affecting the range of interpolation are 

• Curvature of scan lines in the (x , y) system 

• Parallelism of scan lines over the interpolation range 

The first effect was investigated by constructing great circle arcs of varying 
length and orientation, mapping points on these arcs into the data base, and 
finally comparing the mapped positions with positions obtained by interpolation 
between the ends of the arc. Interpolated positions were required to map onto 
calculated positions to an accuracy of 0. 25 array element. This procedure was 
carried out for sets of arcs with origins at the intersections of a 5 x 5 grid 
covering one quadrant. Table 4-2 presents the results of this procedure. 

Non-parallelism of scan lines would restrict the number of scan lines over 
which interpolation ranges. The primary contributor to non-parallelism would 
be the use of the return- scan data in an oscillating mirror scanner system. 

This problem does not exist for coarse resolution data since coarse resolution 
is averaged over a number of scans. For fine resolution, this problem can 
be eliminated by separating the scan lines into odd and even scans and inter- 
polating separately. 
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Table 4-2. Maximum Interpolation Range as a Function of Position 
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To maximize computational efficiency, it is desirable to interpolate over the 
longest range possible. Two schemes for determining this range present 
themselves: 

• Store precalculated range in a table 

• Attempt interpolation to a bench point and reduce range if inter- 
polation is inaccurate 

The precalculation method has an advantage for a system of constant resolution. 
The attempted interpolation method would be simpler to implement and is also 
applicable for variable resolution. 

The largest feasible interpolation region makes a sensible logical block over 
which to perform computations. Such N x M rectangular blocks of data will 
occasionally cross a face boundary. Clearly, interpolation between points in 
differing (x , y) systems is meaningless. Thus a scheme for crossing these 
boundary regions is necessary. 

A simple logical check will identify blocks which cross face boundaries. The 
face number bits in the four-corner benchmarks serial position can be checked 
against each other. If all four corners lie in one face, all points in the block 
must lie in that face; if any corner lies on a differing face, special measures 
must be taken. 

Assuming that the interpolation blocks are small compared to a face, two types 
of face crossing may occur. These are edge crossing and vertex crossing, as 
illustrated in Figure 4-6. 

Two types of vertex crossing may occur (Figure 4-7). Type 1 .is characterized 
by a block having corners on three faces, and this situation may be distinguished 
immediately. Type 2 is characterized by a block having corners on only two 
faces and is difficult to distinguish from a simple edge crossing. To determine 
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whether a block contains an edge crossing or a type 2 vertex crossing, the 
following procedure may be used: 

1. Choose one face coordinate system and calculate the position of all 
four corners in this coordinate system. (Some of these points will 
lie on the extension of the face. ) 

2. Test this quadrilateral to see if it contains a vertex. Two tests 
may be used: 

a. Construct the lines connecting the four corner benchmarks 
and, proceeding around the quadrilateral, determine whether 
or not the vertex lies to the same side of the lines connecting 
the benchmarks. If it does, the vertex is within the block. 

b. A more elegant scheme is based on the fact that all of the 
interpolation blocks will be very nearly rhomboidal. Given 
this, it is then determined whether a point lies within or 
without the rhombus by comparing the areas of the triangles 
formed when the vertex point is joined to two points on the 
rhombus. If any of the four possible triangles has an area 
greater than half the rhombus area, the vertex lies outside 
as illustrated in Figure 4-8. 

When it has been determined which and how many faces of the data base cube 
are involved, the interpolation method then proceeds to calculate the data base 
coordinate system displacement vectors for each face involved. The procedure 
for interpolating across edges (two face) is as follows: (1) starting at the end 
of a scan line, project the end point in both coordinate systems; (2) test to see 
which coordinate system produces a point lying in the face of that coordinate 
system (i. e. , test x and y coordinates against the range of x and y) ; and 
(3) interpolate in that system, testing against the range of (x , y) accordingly. 
When a point exceeds the range of (x , y) , change systems and continue. In 
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the case of two faces, it is not necessary to continue testing the range since a 
scan line will cross only one boundary. The three-face case is similar except 
that, when the test fails, the other two face systems must be tested to determine 
which to use, and testing aginst the (x , y) range must continue until two-face 
boundaries have been crossed. 
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5.1 SUMMARY 


In considering the design of any large-scale data base system, the efficient 
design of I/O management facilities is of great importance. This is especially 
critical in the case of a data base of this nature and magnitude. Up to this 
point, the data base has been considered in a general sense. This section dis- 
cusses a specific data base with properties approximating an operational data 
base. 

5.2 EXEMPLARY MODEL 

Consider a data base containing whole Earth information for one spectral band 
consisting of one data item per pixel. As a proposed model, let each item be 
represented by one byte (6 or 8 bits) and let each scan line contain 1280 pixels 
at a resolution of 1 nautical mile, constant along the scan line. 

The area of the Earth in square nautical miles is approximately 

I 3 \^ 8 

477 X \3.48 X 10 / = 1.52 X 10 

The nearest level of division giving coarser resolution is the 12th level 

6 x 4 12 = 6 X 4096 2 = 100663296 = 1.01 X 10 8 
The linear resolution of this data base is 

= 1. 22 nautical miles 
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8 

Therefore, at this level, 10 bytes of storage will be required. A convenient 
(but not unrealistic) model for disk type storage can be chosen to have six spin- 
dles, 256 cylinders/spindle, 16 tracks/cylinder, and 4096 bytes/track. Thus, 

g 

there will be a total of 6 x 256 x 16 x 4096 = 10 bytes of storage. By splitting 

the serial data string over such a facility, it is advantageous to choose one face 

to correspond to one spindle, and a natural record to correspond to one track. 

Hence, each face would be divided into 256 x 16 = 4096 physical records, and 

6 

each record would represent division of a face at the sixth level (4096 = 4 ). 
Because there are two bits associated with each level of division, a typical 
serial string address would then appear as follows: 
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Each physical record would thus represent an area 64 units square. This size 
is of the same order as the typical propagation distance for the interpolation 
described in Section 4. Figure 5-1 illustrates the relation of various ratios 
of interpolation block size to I/O record size. 

In the case where the two sizes are not identical, and where the blocks and 
records are not perfectly aligned, a certain amount of redundant I/O transfer 
is inevitable. The amount of redundant transfer decreases as the ratio of in- 
terpolation to record size increases. On the other hand, the number of indi- 
vidual I/O operations will rise if this ratio is obtained by lowering the record 
size. The optimum value for this ratio will have to be determined by studying 
the actual machine configuration on which the data base is implemented. 

For example, a 64 x 64 unit record has been chosen to provide reasonable 
savings in computation. Similarly, it is preferable to adopt the same size 
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interpolation block so that there will be 64 x 64 = 4096 interpolations per 
mapping computation. This ratio is illustrated in Figure 5-l(b). To assure 
complete coverage of the interpolation block, six records will be required to 
minimize I/O operations. Consideration must also be given to sufficient core 
buffer space to store all six records simultaneously, requiring 24,576 bytes 
of storage. Data will be queued on a first-in first-out (FIFO) basis, meaning 
that when a record not currently in core is required, it will replace that record 
which has been in core the longest. Figure 5-2 illustrates the typical result 
of such a scheme. Records 1 through 6 are resident in core at one time. When 
record 7 is required, it replaces record 1; record 8 replaces record 2; and 
so on. 

The implementation of such a scheme is accomplished via the use of a push- 
down stack of fixed depth. The stack consists of a set of pointers to the initial 
location of a record's core location and the serial number of the record. When 
a record is required, the stack is searched by record number. If the record 
is present, the associated pointer locates the record; if the record is not pres- 
ent, its serial number is added to the top of the stack. The record which is 
pushed out of the bottom of the stack is returned to disk storage and the new 
record required is read into core in the newly freed locations. The pointer 
which was associated with the oldest record, now is associated with the new 
record. 

The procedure for accessing the data base thus breaks down into two parts: 

(1) the calculation of serial addresses, and (2) the retrieval and storage of 
records implied by the serial addresses. As discussed previously, addresses 
are calculated in blocks whose size is based on the interpolation range. Taking 
a scan line of length 1280 pixels and the chosen 64 x 64 pixel interpolation block, 
it is required to have 64 x 1280 = 81,920 bytes divided into 20 64 x 64 pixel 
blocks. Looping over pixels and scan lines in these blocks, calculations will 
be performed for the serial address as previously discussed. Each serial 
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Order of First Appearance for I/O Records 
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address, as it is computed, is compared under a mask with the previous serial 
address to test for a change in physical record number. If there is no required 
change in record number, the data item is stored into the core location found 
by summing the within-record portion of the address and the base address cur- 
rently in use. On the other hand, if there is a change in the record number, 
the new record desired is compared with the list of records in core (the stack), 
and an I/O operation is performed when appropriate. The base address for this 
new record (as stored in the stack) now becomes the active base address, and 
the data item is stored in the appropriate location. 

The primary features of the preceding scheme are 

1. I/O takes place only when a location is not available in the core 
buffer 

2. All memory references may make use of a base address register 

3. Testing for correctness of the current record requires only a simple 
test under mask operation 

Examination of Figure 5-2 shows that record changes take place approximately 
once per scan line or once per 64 address computations. Thus, each interpo- 
lation block requires access to two or three records not currently in the stack 
(say 2. 5 as a conservative estimate). The time required for a record change 
not necessitating an I/O operation is negligible (six comparison operations 
every 64 address calculations). However, the time required when I/O takes 
place is considerable. 

To calculate the I/O time required, consider the following realistic but con- 
servative set of specifications for a disk storage device: 

• Data transfer rate: 500, 000 bytes/second 

• Average rotational latency: 10 msec 

• Average seek time: 30 msec 
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At the 500K baud transfer rate, one 4096 byte record would take 8 msec. The 
near-neighbor relationship established by the serial addressing scheme assures 
that approximately 75 percent of the immediate neighbors of a record will lie 
on the same cylinder, thus giving an average access time of 15 msec. Using 
the estimate of 2. 5 data base accesses per 64 x 64 block and the fact that each 
access involves one input and one output operation, it is found to require an I/O 
time of 5 x (15 + 8) = 115 msec per block or 2. 3 seconds per 64 scan line input 
block. This yields approximately 900 I/O seconds per orbit, comfortably faster 
than a real-time value of 6000 seconds. 

Another way of looking at this I/O requirement is that it would take approxi- 
mately 2700 seconds or 3/4 hour to fill six faces of 4096 records each. These 
figures are conservative and could be improved considerably by using faster 
hardware and by spreading adjacent records over corresponding cylinders of 
all available spindles. In any case, the I/O requirements are clearly well 
within the state of the art and would not unduly overload a large-scale, multi- 
purpose computer system. 

5.3 HIGHER RESOLUTION DATA BASES AND USE OF SELECTED AREAS 

g 

A whole Earth data base at one mile resolution requires 10 bytes of storage 
per spectral band. This number is within the range of current disk storage 
devices. A higher resolution data base would, however, rapidly exceed the ca- 
pacities of even the largest disk type storage devices. A data base at 0. 1 mile 
resolution would require 10^ bytes, corresponding to literally hundreds of disk 
packs. Mass storage devices, as described in Appendix C, would allow storage 
of such an array. 

The basic record addressing scheme discussed in Section 5.2 requires online 
access to a full six-faced data base. For purposes of storing high resolution 
data, some scheme for storing selected areas is necessary. The simplest 
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method of accomplishing this is by using an indirect addressing scheme. Fig- 
ure 5-3 shows a face divided at the fourth level. 

Blocks of records corresponding to these fourth-level areas can be stored for 

regions of interest. An array containing a logical flag specifying the existence 

of records would be stored along with a pointer indicating the real starting rec- 

2 

prd number for the block. For a fourth-level division, two 16 = 256 arrays 
will be required to store the flags and pointers. 

The I/O management module in the program would first do a shift to determine 
the storage block number, check this against the flag array, and accordingly 
access the record by adding the within-block portion of the record number to 
the starting record number for that block. 
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Figure 5-3. Face Divided at Fourth Level for Windowing 
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SECTION 6 - DATA RETRIEVAL 


6.1 SUMMARY 

This section is concerned with accessing data from the data base for the pur- 
poses of (1) numerical processing and (2) visual display. These are considered 
to be the two main needs of the contemplated user who would be interested in 
using the system. 

6, 2 DATA ACCESS FOR USER NEEDS 

For the two purposes mentioned above, retrieval in physical record units would 
usually make sense. For most visual purposes, shape distortion would be seri- 
ous only in the vertex regions where 120-degree angles map into 90-degree 
angles. This distortion could be minimized by using an overlay grid. 

For numeric purposes, the equal-area property of the data base should make 
the face coordinate systems extraordinarily useful. Most applications will, 
however, require a conversion to a more normal two-dimensional mapping 
(raster) scheme. An algorithm which accomplishes this conversion is pre- 
sented in the source listing attached as Appendix E.9. 

For visual purposes, a useful display scheme would be to display an unfolded 
cube as in Figure 1-6, featuring continental outlines and a data base coordinate 
grid. The user could access any record or group of records by indicating the 
corners of the desired area with a light pen or some other device. The system 
would then retrieve the indicated records, convert to raster form, and display. 
Precalculated geographic or other grids could be overlayed on the data to pro- 
vide referencing. 
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Alternatively, the data base system could be reconverted into a geographic grid 
system by reversing the mapping and interpolation procedure and applying the 
transformation: 


4 = f(x, y) 

V = f(y, x) 

Interpolation and I/O considerations would be similar to those described pre- 
viously except that no output operations would be necessary. 

A particularly rapid scheme for display would involve the use of a special piece 
of hardware, similar to the x, y register summing device discussed earlier. 

A serial address decoder accessing storage via a direct memory access (DMA) 
type device would allow display of a selected set of records with no central 
processor activity other than I/O. In a distributed system, such as that de- 
scribed in Appendix D, this device would be part of a stand-alone display device 
accessing data on a staging device. 
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SECTION 7 - RESULTS OF STUDY 


7.1 SUMMARY 

The primary result of this study is the conclusion that implementation of a 
Quadrilateralized Spherical Cube Data Base is feasible using current large- 
scale computer systems. Numerical results supporting this ponclusion are 
presented below. 

7.2 NUMERIC RESULTS OF MAPPING FUNCTION 

The transformation function f(x , y) described in Section 2.2 was studied in- 
tensively. A FORTRAN program was written to minimize the residual function 
0 . The function was evaluated at 36 locations on a regular grid covering one 
quadrant of a face. (The symmetry of Equation (2-14) permits the generali- 
zation to the other three quadrants of the face.) 

Table 7-1 illustrates the very nearly equal-area nature of the projection. The 
table presents the ratio 


(M _ M *H.) 

\Qx dy dy dx I 



which appears in Equation (2-26). In other words, it presents the transformed 
area at a given point normalized to the projected area at the center of a face. 
For an exact transformation, this ratio is identically equal to unity. The cal- 
culated root mean-square deviation over the entire face is found to be 0.000103. 

The function f*(£ , r\) was obtained from (f(x , y) as indicated in Section 2.3. 
The procedure employed was to minimize the residual function 0* calculated 
on the same regular grid covering one equadrant of a face. 
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Table 7-1. Tabulation of Direct Mapping Function Accuracy 


CENTER 



0.0 

0.2 

0.4 

0.6 

0.8 

1.0 

0.0 

1.0 

1.0001 

0.9999 

0.9999 

1.0001 

1.0 

0.2 

1.0001 

1 .0002 

1.0001 

1.0001 

1.0001 

1.0002 

0.4 

0.0 

1.0001 

1.0 

1.0001 

0.9999 

1 .0002 

0.6 

0.9999 

1.0001 

1.0001 

1.0 

0.9999 

1.0001 

0.8 

1.0001 

1.0001 

0.9999 

0.9999 

1.0001 

0.9998 

1.0 

1.0 

1 .0002 

1.0002 

1.0001 

0.9998 

1.0 


VERTEX 



Table 7-2 presents the results of this minimization, showing the values of 


/([x - f*(f(x, y), f(y, x))] 2 + [y - f*(f(y, x), f(x, y))J 2 j 

over the grid. For an exact inverse transformation, these values are identi- 
cally equal to zero. The RMS deviation of these values is found to be 0. 012. 

7.3 NUMERIC RESULTS FROM TRANSFORMATION AND FAST-FILLING 
ALGORITHM TESTS 

A FORTRAN-coded version of the {£ , V) ( x . y) transformation algorithm 
was tested. This algorithm requires the following arithmetic operations to map 


one point: 

Arithmetic comparisons 3 

Integer multiplication 1 

Integer subtraction 1 

Real multiplications 51 

Real divisions 5 

Real additions 26 

Real subtractions 2 


In this form, the algorithm requires approximately 150 Msec to execute on the 
IBM S/360-91 used for these tests. 

The interpolation and tests used in the simple fast- filling method require 
30 m sec per point to calculate the serial address and store one point into a data 
base record. This time could be reduced somewhat by coding in assembler 
language, but is reasonably close to the irreducible minimum. 

The two vital figures in this study are the 30 Msec calculation time per data 
point and the average I/O time per point. Using the value of 115 msec per 
64 x 64 block obtained in Section 5, this is also approximately 30 Msec per 
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Table 7-2. Tabulation of Inverse Mapping Function Accuracy 


CENTER 



0.0 

0.2 

0.4 

0.6 

0.8 

1.0 

0.0 

0.0 

0.003 

0.007 

0.020 

0.008 

0.0 

0.2 

0.003 

0.008 

0.009 

0.019 

0.016 

0.027 

0.4 

0.007 

0.009 

0.016 

0.011 

0.015 

0.005 

0.6 

0.020 

0.019 

0.011 

0.013 

0.015 

0.025 

0.8 

0.008 

0.016 

0.015 

0.015 

0.008 

0.018 

1.0 

0.0 

0.027 

0.005 

0.025 

0.018 

0.0 


VERTEX 
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point. Overhead from other operations should not contribute significantly since 

3 

all operations will be done in blocks of the order 10 . For a 1280 point scan 
line, computations and I/O times of the order of 0. 02 second per scan line are 
expected. 

For a 100-minute period satellite, the ground track speed is about 3.5 nautical 
miles/second. At 1 mile resolution, this means computation would be taking 
place at about 15 times the data rate. At 0. 1 mile resolution, computation 
would be 1. 5 times faster than real time. 

The above figures are approximate and refer to execution times on the IBM 
S/360-91; they are conservative and could be improved upon by utilizing special 
hardware or by maximizing coding efficiency. They indicate that the imple- 
mentation of a Quadrilateralized Spherical Cube Data Base is feasible on current 
generation, large-scale computers, particularly at coarser resolutions. 
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APPENDIX A - DETAILS OF MAPPING FUNCTION 


A. 1 DETAILS OF DIRECT MAPPING FUNCTION 


Without any loss of generality, Equation (2-22) may be re-written as 


£ = f( x » 


y) =yx + 


(1-7) 3 
x 

r 

o 



2 

y 


g( x , y) 


2 

+ x h(x) 


(A-l) 


where g(x , y) and h(x) are even functions of their arguments. Differentia- 


tion with respect to x yields 


a L r , 3(1 -n 

Sx 2 

r 

o 


2 2/2 
x + y (r Q 



g(x, y) 


2 2 

- 2x y g(x, 


y) 


2 / 2 
+ x y \ r o 


2 

x 


3 g 2 

O X 




h(x) 


4 

2 x h(x) 


(A -2) 


3 

+ x 




dh 

dx 


Substitution of the condition in Equation (2-19) into the above equation yields 


M = 


H 

3x 


x =r c 

y=r c 


3 - 2y - 2 r 


4 

o 


« (r o- 


4 

r ) - 2 r h(r ) 
o o o 


(A -3) 


For convenience, let 6 and to be defined by 


6 - g(r Q , r o ) 


(A-4) 


to = h(r ) 
o 


A-l 



It follows that g(x, y) and h(x) can be expressed as 


g(x, y) 


6 + 




p(x, y) + 




S(x, y) 


(A-5) 


h(x) = LO + 



(A -6) 


where p(x, y) , q(x) , and s(x, y) are even functions of their arguments. 
Moreover, from Equations (A -3) and (A-4), CO can be expressed in terms of 
6 by the equation 





^3 - 2y 


- M - 



(A-7) 


Next, at the point (£ = x = r Q , V = y = 0) , it is noted that Equation (2-12) 
yields 



9 77 

SC 

. *1L 

Sx 

9y 

" Sy 

9x 


x - r 0 

x ~ r o 

x ~ r o 


y=o 

y=o 

y=o 



(A -8) 


However, from the symmetry and boundary conditions imposed on the mapping 
function f(x , y) , it is evident from Figure 1-4 that 



Moreover, from the symmetry of Equation (2-14) 


d rj 

Sy 


x=r 


o 


y=° 


a_£ 

Sx 


x=0 

y =r o 


(A-10) 


Consequently, Equation (A-8) becomes 


M 

a x 

X=r o 

y~o 


a x 


x=0 



y=r 


o 


(A-ll) 


Equations (A-2), (A-5), (A-6), and (A-7) yield 


H 

dx 


x~ r c 

y=o 


= 3 - 2y - 2 r h(r ) 
o o 


= 3 -27 - 2r CO 
o 


4 

= U + 2 r 

o 


6 


(A-12) 


a 

a x 


x=0 


4 

= y + r 


o 


g(°. r Q ) 


y =r o 


= y + r 


4 

o 



s(0, r Q ) 


(A-13) 


A -3 



Next, it is noted that the number of terms in g(x, y) given by Equation (A-5) 
will be minimized by choosing 


s(x, y) = 0 


(A-14) 


This will also eliminate the terms s(0 , r) in Equation (A-13), which then 
becomes 


H 

bx 


x=0 


= y + r 6 
o 


y= r o 


(A-15) 


Substitution of Equations (A-12) and (A-15) into Equation (A-ll) yields 


H + 2 r 4 + r 4 6) = 2y/2y* 


(A-16) 


or, equivalently, 


2r 8 6 2 + m + 2y)r 4 6 + yiy - 2\/2y 2 ) = 0 
o o 


(A-17) 


Solution of this equation for r^ 6 yields 


4- -B ± y ( B 2 - 4AC) 

r* 0 = 

o 2A 


(A-18) 


A-4 



where 


A = 2 

B = (fJ, + 2y) > 0 (A-19) 

c = {[iy - 2, fly 2 ) < 0 


Consideration of Figure 1-4 reveals that 6> 0 . Hence, only the plus (+) sign 
is retained in Equation (A-17). Thus, the following expression is obtained 
for 6 : 


6 = 


4 r 


- (/i + 2 y) + J{\? - ±\iy + 4y 2 + Vofiy 2 ) 


= 0. 79048 64491 208 


(A-20) 


The value of CO is obtained by using Equation (A-7): 


1 / „ 4, 

CO = — (3-2y-/i-2r 0 

2 r 

o 


= -1.2254 41487 984 


(A-21) 


Finally, in view of the above discussion, the direct mapping function f(x , y) 
can be chosen to have the form 


„ (1 - y) 3 2/2 2 

f(x, y) = yx + — x + xy r - x 


6 + (r 2 - y 2 ) c. . 


2i 2j 

x y 


i^O 

j-0 


(A -22) 


+ X 3 (r 3 - X 2 ) a> + (r 2 - x 2 ) Yj d i 

L i^O 


2i 


x 
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A. 2 DETAILS OF INVERSE MAPPING FUNCTION 


From Equation (2-12) and the following property of Jacobians 



it follows that the inverse transformation 


x = x(£, V) 
y = y(€, V) 


(A -2 3) 


(A -24) 


satisfies the following equation 


Sx dy 

ac srj 


Sx a y 

a^ a~c 


i 



2\3/2 


(A -2 5) 


Because the direct transformation in Equation (2-14) is symmetrical, the 
inverse transformation is also symmetrical, i.e., 


x = f*(|, V) 
y = f * (r?, 0 


(A -26) 


Proceeding as in Section 2.2, it is seen that the following four conditions are 
the counterparts of Equations (2-16) through (2-19) 

f*(0, V) = 0 (A -27) 


A -6 



(A -2 8) 


f*(r , 7?) = r 
o o 


dx 

H 


C=o 

77=0 


Zy 

bn 


£=o 

77=0 


_ i = y * 

y 


(A -29) 


dx 


£ =r o 


dr? 

Sy 


i= r o 


1 - * 

= — — u * 

JU H 


(A -30) 


7?=r 0 r?=r 0 


Thus, the first three conditions above may be incorporated into the following 
expansion which is the analog of Equation (A — 1 ) : 


x = f*(4, n)=yH + { ^P t 3 + t(v 2 0 - 4 2 


v 2 g*(Z, v) + 4 2 h*(0 


(A -31) 


where 


g*(4, V) = 5* + (r 2 - 7] 2 ] p*(£, *?) + (r 2 - £ 2 ) s*(4, 


71) 


(A-32) 


/ 2 2\ 

h*(£) = to* + ^r - £ j q*(£) 


(A-33) 


6* = g*(r , r ) 
& ' o o' 


to* = h*(r Q ) 


(A-34) 


p*(£, V) , q*(£) » and s*(4, V) are even functions of their arguments. 


A-7 



The fourth condition, i.e. , Equation (A-30), yields 


to* 


—L. (s _ 2y * - 11 * - 2 r 4 6*) 
2 r 

o 


(A-35) 


which is the analog of Equation (A-7). 

Next, by differentiating Equation (A-31) and using Equations (A-32) through 
(A-35), the following analogs of Equations (A-12) and (A-15) are obtained: 


Sx 

H 


4=r 0 

7)=0 


3 


- 2y* - 2 r 


4 

o 


CO* 


(A -36) 


= ji* + 2r 


4 

o 


6 * 


4 

= y * + r 

2 

6 * + r s*(0, r ) 

vf*Y 

II 

O 

c 

o o 


7)=r 

o 


(A-37) 


However, consideration of Figure 1-4 reveals that the following expressions 
are also valid: 


dx 




a 4 


Sx 



4= r o 


o 

f-i 

ii 

X 


7)=0 


y=o J 



+ 2r 



Sx 




-i i 


4=0 

Sx 

x=0 

4 , 

y + r o 
o 


r7=r 0 


y= r o_ 



= y* 


(A-38) 


(A-39) 
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By equating the corresponding expressions above and then solving for 6* and 

s*(0 , r ) , it is found that 
o 


(M* " M*) 

6* = — — - — (A -40) 

2 r 

o 


s*(0, r ) =~ 

O 2 

r 

o 

For convenience, let 6* be defined by 

6* = s*(0, r ) (A -42) 

1 o 

Then, it follows that s*(£, V) can be expressed as 


O'* - y*) 


6 * 


(A-41) 


-V V) (A -43) 

where u*(£) and v*(£ , T}) are even functions of their arguments. For the 
purpose of expanding x = f*(£ , r?) and also satisfying the six conditions de- 
scribed by Equations (A -27) through (A -30) and (A-38) through (A -39), it may 
be verified that it suffices to let 


s*(4, 77) = 6* + i u*(£) + r 


u*(£) = 0 


v*(£, V) - o 


(A -44) 


A-9 



Finally, in view of the above discussion, the inverse mapping function f*(£, r\) 
can be chosen to have the form 


f*(l, V) I 3 + £r? 2 (r 2 - £ 2 )6* + 6*|r 2 - 


/ 2 2\ V - * . J2.i 2jl .3/ 2 J2\ 

+ ( r 0 - , »)L c fj« * + M r 0 - 4 ) 


i^O 

1*0 


CO* 



(A-45) 
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APPENDIX B - DEPENDENCE OF SCAN ELEMENT LOCATION ON 
ORBIT AND ATTITUDE OF SATELLITE 

B. 1 THE CASE OF ZERO PITCH, ROLL, AND YAW 

A A 

Consider an Earth-centered "auxiliary" orbital coordinate system (x , y , 

A 

z ) , as illustrated in Figure B-l, such that the satellite is always located on 

A A A 

the y -axis and the scan angle a is always in the (y , z )-plane. Then, the 

A A A 

z -axis is always perpendicular to the plane of the orbit, and the x - and y - 

A 

axes move within the orbital plane, rotating about the z -axis. 

It is noted that the point of intersection, Q , is such that |p| <77/2 so that 
cos fj, is always given by 


cos /i =J(1 - sin 2 /i) 

From the sine formula for AOPQ , it follows that 


(B-l) 


sin cr sin (77 - a - \i) 
R ” (R + h) 


Using the relation 


sin (77 - ct - /i) = sin (a + /i) = sin cr cos ^ + cos cr sin ^ (B-3) 


Equation (B-2) becomes 


sin a sin cr cos /I •+ cos cr sin {J, 
R (R + h) 


(B-4) 


B-l 






Figure B-l. Auxiliary Orbital Coordinate System 


B-2 


or, equivalently, 


cos H + cot a sin JU = 


(R + h) 
R 


(B-5) 


Using Equation (B-l), this becomes 



2 /D 

sin J M) + cot a sin jj,= — 


h) 


R 


(B-6) 


However, it is noted from AOQS that 

A 

S injj, = £-— (B-7) 


Hence, Equation (B-6) becomes 



or, on squaring, 


or 


R R R 


(B— 9) 


2 z 
(1 + cot a) — 


R 



A 

’ , 2 

„ (R + h) 

2 - — — — - cot c 

z 

_ + 

(R + h) x 

R 

R 

2 

R 


= 0 


(B-10) 


B-3 



Solution of this equation yields 


z A _ -B ± y/(B 2 - 4 AC) 
R 2 A 


(B-ll) 


where 


2 2 

A = 1 + cot a = esc a 


B = -2^-^cot a 
R 


(B-12) 


C - - ! 

L ■ 2 
R 


Since the point of intersection, Q , is on the closer side to P , its z -coordinate 
is given by Equation (B-ll) in which only the minus (-) sign is retained. Thus, 
Equation (B-ll) becomes 


A 

z 

R 


(R + h) 
R 


. 2 

sin o cos o - sm 


2 (R + h) 
esc a - 


R 


(B-13) 


Also, from Figure B-l it is noted that 


A /~2 

y = R cos I J = yj { R 



and 


A 

x 



(B-14) 


(B-15) 


B-4 



Next, consider the Earth -centered inertial coordinate system (x , y , z ) , 
as shown in Figure B-2. 

Then, it is well-known that the Euler angles 0 , 0 , and 0 are given by 

0 = O (B-16) 


0 = i 


(B-17) 


ip = 03 + v - TT/2 (B-18) 

where i = inclination of satellite orbit 
= longitude of ascending node 
03 = argument of perigee 
v = true anomaly 

The coordinates of point Q are denoted by (x 1 , y 1 , z 1 ) in the two coordinate 
systems. For convenience, let 


-oJ - , I I I. 
r - (x , y , z ) 


_^A _ A A A 
r = (x , y , z ) 


I A 

Then, the vectors "r" and *r“ are related by 


(B-19) 


— 1 a -!^ a 

r = A r 


(B-20) 


-1 


where A is the matrix given by Equation (4-47) of Reference 2: 


1 cos ip cos 0 - cos 0 sin 0 sin lb -sin Ip cos 0 - cos 0 sin 0 cos Ip sin 0 sin 0 

A -1 = A = | cos ip sin 0 + cos 0 cos 0 sin tp -sin |() sin 0 + cos 0 cos 0 cos tp -sin 0 cos 0 

sin 0 sin tp sin 0 cos tp cos 0 


(B-21) 


B-5 




Figure B-2. Inertial Coordinate System 


B-6 



Hence, from Equations (B-15), (B-20), and (B-21), it follows that 

I A A 

x = (-sin 0 cos 0 - cos 6 sin 0 cos 0) y + (sin 0 sin 0) z (B-22) 

I A A 

y = (-sin 0 sin 0 + cos 0 cos 0 cos 0) y + (-sin 0 cos 0) z (B-23) 

T A A 

z = (sin 0 cos 0) y + (cos 0 ) z (B-24) 

A A 

where y , z , 0 , 0 , and 0 are given by Equations (B-14), (B-13), (B-16), 
(B-17), and (B-18), respectively, Hence, the coordinates (x 1 , y 1 , z 1 ) of the 
point Q are obtained. 

Finally, consider the uniformly, Earth-centered rotating coordinate system 

73 73 73 

(x , y , z ) as shown in Figure B-3. 

It is easily shown that 


R I 

X = X 


cos Wt + 


I 

y 


sin tot 


(B-25) 


R 

y 


= -x sin tot + y cos tot 


(B-26) 



(B-27) 


in which t = 0 when the two coordinate systems coincide. Thus, the coordi- 
nates of a scan element on a rotating Earth are obtained in terms of the satel- 
lite orbital position, satellite attitude, and scan angle. 
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Figure 33-3. Uniformly Rotating Coordinate System 



B. 2 THE CASE OF GENERAL PITCH, ROLL, AND YAW 

AAA III R R R .. . , 

Let (x , y , z ) , (x , y , z ) , and (x , y , z ) , respectively denote 

the "auxiliary" frame, the Earth-centered inertial frame, and the Earth- 

centered rotating frame of reference as defined previously. 

Consider an Earth-centered orbital coordinate system (x° , y° , z°) defined 
by the unit vectors 


y\0 /\A 

X = -x 



(B-28) 


Figure B-4 illustrates the orientation of the orbital coordinate system with 
respect to the "auxiliary" orbital coordinate system. 


Let 


_A _ A A A, 
r = (x , y , z ) 

(B-29) 

^o _ o o o v 
r = (x , y , z ) 


Hence, from Figure B-4, it is seen that 


r A" 


- 


* o" 

X 


-10 0 


X 

A 




o 

y 


0 0-1 


y 

A 




0 

_z 


. 0 -1 o _ 


_ z 


(B-30) 
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Therefore, the transformation 



has the matrix B given by 



0 

0 

-1 



(B-31) 


(B-32) 


It is noted that the satellite pitch, roll, and yaw axes are obtained by a parallel 
translation of the (x° , y° , z°) coordinate system from O to P . 

B B B 

Next, consider a body-fixed coordinate system (x , y , z ) on the satellite. 
Let a , (3 , and y respectively denote the roll, pitch, and yaw angles which 
transform a vector from the orbital to the body-fixed frame. That is, if 


then 


_^o _ , o o 
r = (x , y , 


z°) 



B 

y > 



(B-33) 



(B-34) 


B-ll 



where 


cos 8 cos y 
+ sin P sin ol sin y 


cos ol sin y 


-sin P cos y 
+ cos P sin cl sin y 


C = 


-cos P sin y 
+ sin 8 sin a cos y 


cos ol cos y sin P sin y 

+ cos P sin a cos y 


(B-35) 


sin j8 cos a 


-sin ot 


cos P cos OL 


Since C is an orthogonal matrix, 



(B-36) 


Hence, from Equation (B-35), the following matrix is obtained for 




Vr 

-c a s 

8 y 

s„c 

p OL 

+ Wr 

+ S 0 S C 
(3 ol y 


c s 
a y 

c c 

ol y 

-s 

OL 

-s R c 
p y 

S_S 

£ y 

C Q C 

P OL 

+C Q s s 
pay 

+c Q s c 

P ol y 



(B-37) 


where the abbreviation S and C are used to denote sin ol and cos ol , 

OL OL 

respectively. 

Let n denote the unit vector along the scan ray which is assumed to be in the 
B B 

(y , z ) -plane. The scan angle, ct , is defined as in Figure B-5. 

Then, in the body-fixed system, the scan ray unit vector is given by 


/\ aB . aB 
n = n = -sin cry 


A B 

+ cos cr z 


(B-38) 
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Figure B-5. Satellite Scanning Earth 
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In general, if the scan ray is not in the (y , z ) -plane, then n is given 
by 


n = n B = cos a x B + cos cr y B + cos cr z B (B-39) 

x y Z 

where cr O' and o are the angles between n and the respective body- 
x ’ y ’ z 

fixed coordinate axes. 

It follows from Equations (B-31) and (B-34), that in the orbital and the auxiliary 
orbital systems, the scan ray unit vector is given by 

= C _1 ft B (B-40) 

= Bft° (B-41) 

so that 

n A = B C _1 n B (B-42) 

From Equations (B-32), (B-37), and (B-38), it follows that 


-c, c 
p y 

C a S 

0 y 

-s a c 

p a 

- s 0 s s 
/3 a y 

- s 0 s c 
j3 a y 


So c 
i3 y 

-S a S 

P y 

-c a c 

P a 

- C,S S 
pay 

- c a s c 
Pay 


-C s 
a y 

-c c 
a y 

S 

a 
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-1 aB 

BC n = 


C a )<C o )} 

U-^ y -c ffSoi C y) (- Su ) + (-c^c a HC a )} 
U-c a c y H-s a ) + V c pc a >) 


(B-44) 


Hence, if we denote n by 


A a A A^A A^A 

n = n x +ny + n z 
x y z 


(B-45) 


it follows from Equations (B-31), (B-33), and (B-34) that 


"x = < c 3 s r-^ s a C r )< 'V + ( - s 3 C a )(C a ) 




(B-46) 


V(- C a v ( ^ + W 


It is obvious from Figure B-5 that by definition 


n A < 0 

y 


(B-47) 


Next, let (x A , y A , z A ) denote the coordinates of the point Q in the auxiliary 
orbital system. It is noted that the coordinates of P are (O , R + h , O) . 
Hence, it follows that the following ratios are equal 


x A y A - (R + h) z A 


(B-48) 


B-15 



which comprise a system of two independent equations and three unknowns. The 
third equation is obtained by noting that the point Q lies on the sphere 



Hence, we obtain 


A 

n 



n 

z 


A 

A y A 

y = -A- z +(R+h) 

n 

z 



A 2 

z + 


n 


V A 
A- z + (R + h) 


n 


+ z 


A 



or, equivalently, 



o (R +h) 

— 

A 

n 

y 

A 

z 

, 2 

(R + h) 

2 R 

A 

~R~ + 

2 ' 1 


n 

z 


R 


(B-49) 


(B-50) 


(B-51) 


0 (B-52) 
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It is easily verified that Equation (B-52) reduces to Equation (B-7) for the case 
of zero pitch, roll, and yaw. Solution of Equation (B-52) yields 


z^_ -B ±J( B 2 - 4 AC) 

R 2 A 


(B-53) 



(B-54 


Since the point of intersection Q is on the closer side to P , its z -coordinate 

is given by Equation (B-53) in which only the minus (-) sign is retained. Equa- 

A A 

tions (B-50) and (B-41) are then used to yield x and y . 

The rest of the analysis is precisely the same as that in Section B. 1 commenc- 
ing with Figure B-2. Thus, the coordinates (x* , y 1 , z 1 ) of Q in the inertial 
R R R 

system and (x , y , z ) in the rotating system are readily obtained by the 
use of Equations (B-22) through (B-24) and (B-25) through (B-2 7). 
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APPENDIX C - USE OF MASS STORAGE DEVICES 
AND THREE-DIMENSIONAL ARRAYS 


The discussion of data bases in this report has assumed that only current data 
would be maintained in direct-access storage. A new class of storage devices, 
rapid mass stores, are now becoming available. These devices open up the 
possibility for on-line access to noncurrent data. 

These devices are characterized by the ability to store on-line up to trillions 
12 

(10 ) of bits of information, and to transfer large blocks of data at high rates 

exceeding 10 bits/seconds. A wide variety of mechanical and electro/optical 
devices with these characteristics is presently coming onto the market. These 
devices offer either write-once/multiple-read or multiple -write and read 
capabilities. 

A write-once system would be suitable for archiving meteorological data bases 
for later reference, and for processing of temporal variations. Data can be 
transferred to the write-once device in the case of data bases with a random 
access device. There are four ways of archiving the data: 

1. Storing raw data for later insertion into data base format 

2. Storing a snap dump of the data base on completion of a satellite 
orbit 

3. Storing a snap dump on a sparse schedule 

4. Storing data in a three-dimensional data base (including a temporal 
coordinate) 

This last method may have excellent long-range utility in that data would be 
blocked together in the temporal as well as spatial sense. This permits rapid 
access to temporally varying data in small spatial regions without long-range 
searches in the mass store. This is important since all mass storage devices 
utilize some sort of mechanical transport for units of the storage medium. 


C-l 



A three-dimensional data base can be constructed by extending the serial 

stringing scheme to 3 bits per level, utilizing 1 bit for each coordinate. A 

non-cubic configuration (with fewer temporal than spatial intervals) can be 

achieved by adding the third bit only to N levels of division, the number of 

N 

temporal intervals being 2 . The natural temporal unit would be one orbit 
which is approximately 100 minutes. The third bit could be added either at the 
highest or lowest level, depending on the application, for efficient access to 
areas selected by temporal or spatial order. 

12 

A device with 10 bits would have the capacity for storage of more than 

g 

100 complete 10 byte data bases, and additional intervals for usage on se- 
lected spatial areas. 

One mass storage device familiar to CSC personnel is the UNICON device used 

12 

in the ILLIAC IV system. This device gives on-line access to 2. 9 x io bits. 
The storage medium is mylar strips which are written by burning areas with 
a focused laser beam. This device thus belongs to the write-once class. 

Mechanically, the UNICON resembles the IBM 2301 Data Cell device, utilizing 
a rotating read/write drum on which the storage strip is positioned. Automatic 

mechanical pickers give on-line access to 1000 storage strips, each containing 

9 4 

2. 9 x 10 bits of data. The data are spread over approximately 10 tracks/strip 

5 6 

with a density of about 2.5 x 10 bits/track. Data can be transferred in 10 -bit 

0 

blocks at a rate of 4 X 10 bits/second. The interblock access time is 
200 msec. 

The general I/O management scheme described in Section 5 would still apply 
to such a device but data could be transferred in much longer records. 
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APPENDIX D - USE OF DISTRIBUTED PROCESSING 


Implicit in Section 6 was the idea that the data base system would be imple- 
mented on a classic large-scale computer system. Substantial savings in 
central processor usage and greater flexibility could be obtained by distributing 
a number of system functions into various subsidiary devices. 

A straight large-scale computer implementation would appear as shown in 
Figure D-l. 

A simple improvement would be to place all I/O management under the control 
of a minicomputer which would handle all data base accesses and transfer data 
to the CPU via a DMA, extended core, or a staging drum as illustrated in 
Figure D-2. 

A further elaboration would involve more elaborate minicomputer systems 
which would handle data input and user interaction functions. Figure D-3 shows 
an arrangement where a system of minicomputers relieves the CPU of all oper- 
ations except calculations based on data transmitted from the data base. In 
such a system, either the CPU or the control minicomputer could control the 
system. Functions indicated as being performed on one minicomputer (such 
as input processing and mapping to the data base system) could actually be 
spread among a number of processors. 

In the input processor area, orbit determination and attitude determination 
functions could take place in parallel in two or more units. In the mapping 
process, the following areas are amenable to parallel calculation: 

• Mapping the £ , rj systems to x , y may be done in parallel 
£ -» x and 77 y 

• Any number of interpolation blocks or scan lines could be processed 
in parallel if serial coordinates were outputted as strings to a 
separate I/O processor 
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I/O processing for any number of spectral bands could be carried 
on in parallel 
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Large-Scale Computer Implementation With I/O 
Management Under Minicomputer Control 
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Figure D-3. Large-Scale Computer Implementation Where 
Either CPU or Minicomputer Controls System 








APPENDIX E - LISTINGS OF TEST AND EVALUATION PROGRAMS 


E.l SUMMARY 

This appendix contains listings of FORTRAN programs and subroutines used 
for testing and evaluating the various transformations and procedures covered 


in this report. 

Listing Section 

Program used for evaluating coefficients of f(x , y) E.2 

Program used for evaluating coefficients of f*(£ , 77 ) E.3 

Mi nim i z nr subroutine used by the programs of Sec- E.4 

tions E.2 and E.3 

Subroutine used for generating a great circle arc E. 5 

Subroutines used for mapping a point in Cartesian E. 6 

space into data base coordinates 

Subroutines used for simulating the fast-filling method E. 7 

Subroutines used for testing I/O handling method E. 8 

Subroutines used for retrieving data from the data base E. 9 

and for graphic display 
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E. 2 PROGRAM USED FOR EVALUATING COEFFICIENTS OF f(x , y) 


1001 


101 


MAIM 

THIS PROGRAM COMPUTES COEFFICIENTS OF THE DIRECT TRANSFORMATION 

FUNCTION RY MINIMIZING THE DIFFERENCE IN AREA FOR SMALL 

ARFAS 

PROGRAMMER - E. MICHAEL O'NEILL 

COMPUTER SCIENCES CORPORATION 

8728 COLESVILLE RD. SILVER SPRING, MD. 

IMPLICIT R E A L * 8 (A-H,0-Z) 

EXTERNAL FUNK 

COMMON RM , R0 , GAMMA , OMGAM,XM,YM, XI N, YIN, DX,DY, DELTA, OMEGA 
COMMON /BILBO/ CHI SO , X ( 20 > , XM I N ( 20 > , XMAX ( 20 ) , DELTAX ( 20 > , 
fi DFLMIM ( 20) ,FRR ( 20 ,20) , MV, NTRACE , MAT R I X , MASK (20) 

COMMON /FRODO/ NFMAX , NFL AT , JV ARY ,NE XTR A 
DATA PI/3. 141592653589D0/ 

DATA NOUT2/10/ 

SFT PRECALCULATED COEFFICIENTS 
DE L TA= 0. 7904 864 A 9 120 8D0 
C,AMS0=PI /6.0D0 
GAMMA=DSORT (GAMSO) 

0MGAM=1 .ODO-GAMMA 

AMU=DSORT (DSORT (3.000) *P 1/ 2.0D0) 

SET DIMENSIONAL PARAMETERS 
RM=1 .0D0 

R0=DS0RT( RM*RM/3.0D0) 

R2= R0*R0 
R4=R2*R 2 

DMFGA=0.5D0/R4* (3 ,0D0-2.0D0*G AMM A-A MU-2. 0D0*D ELT A*R4) 

XM=R0 
YM = R0 

X I N = X M / 5 . 0 DO 

Y I N= X I N 
DX=XM*1 .00-1 2 
D Y = DX 

WRITE! MOOT 2 , 1001 ) 

FORMAT!' INITIAL CONDITION OUTPUT,' 
l'ALL VARIABLF PARAMETERS ZEROED' / ) 

COM PUT F AND DISPLAY INITIAL AREAS 
CALL 0 I SPS 0 ! 6 ) 

CALL D ISPSO (N01JT2 ) 

CONTINUE 

READ INITIAL VALUES AND LIMITS FROM CARDS 
READ! 6, 1 , E ND = 99 ) NV, NFMAX, NTRACE, MATRIX, MASK 
WRITE ( 6,2 ) N V , NFMAX, NTRACE , MATRIX, MASK 
F0RMAT(4I4,4X,20I1 ) 


MATRIX, MASK ' , 4 1 4 , 4X , 20 1 1 ) 


FORMAT! ' 'NV, NFMAX, NTRACE 
READ (5,5) ( X ( I ) , I = 1 ,NV ) 

WRITE (6,6 ) (X ( I ) , I =1 , NV ) 

READ! 5,5 ) ( XMAX! I ) , 1=1 ,NV) 

WRITE (6,6) (XMAX ( I ) , I =1 ,NV) 

READ! 5,5 ) ( XM IN ( I ) , 1 = 1 ,NV ) 

WRITE ( 6,6 ) ( XM IN ( I ) , I =1 ,MV ) 

READ (5,5) (DELTAX! I ) ,1 = 1 ,NV ) 

WR I TE ( 6 ,6 ) (DELTAXII ) , I =1 , NV ) 

RE AD ( 5 , 5 ) (DELMIN! I ) , 1 = 1, NV) 

WRITE (6,6) (DELM IN! I > , 1=1 ,NV) 

5 FORMAT! 10F8.3 ) 

6 FORMAT! 1X,10F8.3) 

CALL DISPS0I6) 

CALL STFPIT TO MINIMIZE RMS AREA DEVIATIONS AS CALCULATED 
BY SUBROUTINE FUNK 

CALL STEP IT ( FUNK ) 

C DISPLAY RESULTS 

CALL 0 I S P S Q ( 6 ) 

WR I T E ( NO) IT 2 , TOO? ) CHI SO , ( X ( I ) , I = 1 ,NV ) „ , 

mo? FORMAT!' CHISn,'/ARIAPLF PARAMFTERS ’, 012,5/(3021.14/)) 
CALL niSPS0(N0UT2) 

GO TO 101 
99 CONTINUE 
STOP 
END 

BLOCK DATA 

C INITIALIZE STEPIT COMMON AREAS 

IMPLICIT REAL*8 (A-H,0-Z) 

COMMON /BILBO/ CH I SO , X ( 20 ) , XM IN ( 20 ) , XM AX ( 20 1 1, D ELT AX ( 20 ) , 
CDELMIN(20! ,ERR(20»20) ,NV, NTRACE, MATRIX, MASK(ZO) 

DATA X/20Y0.0D0/, XMAX/20*0. ID- 10/, XMIN/20Y-0. ID - 10/, 

* DFLTAX/20*0.1D-11/, DELMIN/20*0.1D-13/ 

END . " 


00000300 
00000400 
00000500 
00000600 
00000700 
00000800 
00000900 
0000 1000 

00001100 

00001200 

00001300 

00001400 

00001500 

00001600 

00001700 

00001800 

00001900 

00002000 

00002100 

00002200 

00002300 

00002400 

00002500 

00002600 

00002700 

00002800 

00002900 

00003000 

00003100 

00003200 

00003300 

00003400 

00003500 

00003600 

00003700 

00003800 

00003900 

00004000 

00004100 

00004200 

00004300 

00004400 

00004500 

00004600 

00004700 

00004800 

00004900 


00005000 

00005100 

00005200 

00005300 

00005400 

00005500 

00005600 

00005700 

00005800 

00005900 

00006000 

00006100 

00006200 

00006300 

00006400 

00006500 
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SUBROUT INR DISPSO(NOUT) 

THIS SUBROUTINE PRINTS OUT THE ABSOLUTE AND RELATIVE AREAS 
OF A GRID OF SMALL SQUARES IN ONE OCTANT 

IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION S A ( 6 * 6 ) « AA(6,6) 

DATA SA, AA/36*O.ODO,36*O.ODO/ 

CALL CNTARV(AA) 

SAC=AA( 1,1) 

DO 30 J = 1 » 6 
DO 30 1=1 , J 

30 S A ( I , J ) = AA ( I , J ) /SAC 

WRITE! N OUT ,40) (<AA(I,J),I=1,J),J=1,6> 

WRITE (MOOT, 45) 

WR I T E ( NOIJT » 40 ) (<SA(I,J),I=1,J),J=1,6) 

WRITE ( NOUT ,45 ) 

WRITE! N OUT. AS) 

AO FORMAT! 1X,D12.5/1X,2D12.5/1X,3012.5/1X,4D12.5/1X,5D12.5/ 

1 1 X * 6D1 2 . 5 ) 

A 5 FORMAT! 1H ) 

RETURN 

END 

SUBROUTINE FUNK 

THIS SUBROUTINE CALLS C NT ARY TO CALCULATE THE PROJECTED 
AREAS OF A GRID OF SMALL SQUARES. 

THF RMS OF THE DEVIATIONS OF THE AREAS FROM THE AREA CALCULATED 
AT THE CENTER IS RETURNED AS CHISO 


IMPLICIT REALA8 (A-H,0-Z) 

COMMON RM,RO,GAMMA,OMGAM,XM,YM f XIN,YIN,DX,DY 
COMMON /BILBO/ CH I S 0 , X ! 2 0 ) , XM I N ( 2 0 ) , X MAX ( 20 > , D EL TA X ( 20 ) 
£ DFLMIN (20) .ERR ( 20 ,20) ,NV,NTR ACE , M AT R I X , MA SK ( 20 ) 
DIMENSION SA ( 6 * 6 ) » AA(6,6) 

EQUIVALENCE <C01,X(1)),(C10,X(1)),(C11,X<2)), (C02, X! 3) ) , 
#* ( C12,X (4) ) , (C21 ,X(4> ) , (C22,X<5> ) , < C03,X(6> > * (C30, X<6> ) , 
*, (C13 ,X (7) ) , (C23,X (8) ) , <C32,X(8) ) , (C33»X(9) ! 

CALL CNTARY(SA) 

S AC= S A ( 1 » 1 ) 

DO 20 J=l,6 
DO 20 1=1, J 

20 SA! I , J ) = SA ( I , J ) /SAC 
DEVSO=O.ODO 
DO 30 J= 1 , 6 
DO 30 1 = 1, J 

30 D F V S 0= 0 E V S 0+ ( S A ( I , J ) -1 .000) **? 

CHISO=OSORT (DEVSQ/ 21.0D0) 

RETURN 

END 


00010000 


00010100 

00010200 

00010300 

00010400 

00010500 

00010600 

00010700 

00010800 

00010900 

00011000 

00011100 

00011200 

00011300 

00011400 

00011500 

00011600 

00011700 

00011800 

00016300 


00016400 

00016500 

, 00016600 
00016700 
00016800 

( C20 , X ( 3 ) 100016900 

(C31,X(7) 100017000 
00017100 
00017200 
00017300 
00017400 
00017500 
00017600 
00017700 
00017800 
00017900 
00018000 
00018100 
00018200 
00018300 
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5 

10 


20 


SUBROUTINE CNTARY(SA) 

THIS SUBROUTINE CALCULATES THE PROJECTED AREAS OF A GRID OF 
SMALL SQUARES BY FIRST APPLYING THE DIRECT TRANSFORMATION 
TO A SET OF POINTS, COMPUTING THE AREA OF THE TRANSFORMED SQUARE 
AND PROJECTING TO THE SURFACF OF THE SPHERE 

IMPLICIT REAL*8 (A-H,0-Z) 

COMMON RM, RO, GAMMA, OMGAM,XM,YM,X IN, YIN,DX,DY 
DIMENSION P 1 t 2 ) , P 2 < 2 ) , P3(2), PA ( 2 ) , SA(6,6) 

SPHERICAL PROJECTION FUNCTION , 

SPROJI A,B,C,D,E )=A*( B*B*C) /( ( C*C+D*D+E*E )**1.50D0) 

x=o.ono 
DO 10 J= 1 , 5 
Y=0.000 
DO 5 1=1 
DIRECT 
CALL XY 
CALI. XY TO PF 
CALL XY TO PE 
CALL XY TO PE 
COMPUTE AREA ... 

S1=ARFA(P1,P2,P3) 

S2= AREA (PI ,P4,P3) 

S PRO J FC T TO SPHERICAL SURFACE 
SA( I , J ) = SPR0J(S,RM,R0,P1( 1),P1(2)1 
y = Y + Y IN 
X= X+X IN 
X=0 .ODO 
Y=YM 

DO 20 1=1,6 


f RANSFURMaT IUIm 

TO PF (X,Y , P.o, PI ( 1) ,P1 ( 2) ) 

( v+nx , v ,P0,P2( 1) , P 2 ( 2 ) ) 
(X+DX,Y+DY,R0,P3(1) ,P3(2) ) 
(X , Y + DY ,R0,P4 ( 1 ) ,P4 ( 2) ) 
OF TRANSFORMED SQUARE 


CALL 

CALL 

CALL 

CALL 


XY 

XY 

XY 

XY 


PF 

PF 

PF 

PF 


( X 
( X- 


,Y 

■OX,Y 


$A(I, .. 
X=X+XIN 
CONTINUE 
Return 
EN n 


, R 0 , P 1 ( 1) ,P1(2) ) 
v a — t > a , i , R 0 * P 2 ( 1 ) , P 2 ( 2 j j 

( X-DX , Y-DY , R 0, P3 ( 1 ) , P 3 ( 2 ) ) 
TO PF (X , Y-DY .RO, PA ( 1 ) . P4 ( 2 ) ) 
= ARE A (PI ,P2,P3) + ArEA(P 1,P4.P3) 
ft) = SPROJ(S,RM,RO,Pl (1) ,P1(2! ) 


TO 

TO 

TO 

TO 


00006600 


00006700 

00006800 

00006900 

00007000 

00007100 

00007200 

00007300 

00007400 

00007500 

00007600 

00007700 

00007800 

00007900 

00008000 

00008100 

00008200 

00008300 

00008400 

00008500 

00008600 

00008700 

00008800 

00008900 

00009000 

00009100 

00009200 

00009300 

00009400 

00009500 

00009600 

00009700 

00009800 

00009900 
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SUBROUTINE XY TO PE ( X I , Y I , R , P , E ) 
DIRECT TRANSFORMATION FUNCTIONS 
IMPLICIT RFAL-8 (A-H,0-Z) 


COMMON RM,RO, GAMMA ,OMGAM,XM ,YM ,X I N r Y I N , DX r DY , DFLT A , OMEGA 
. . - - - ,X(20) ,XM IN t 20 > .XMAXI20) ,DELTAX(20) , 

' NTRACE, MATRIX, MASKI20) 


COMMON /BILBO/ CHISO 

£ DELMIN(20),ERR(20.20)»NV}i'iiKm,c,nAiKiAti“iA;>iM 
EQUIVALENCE (C00,X(1)),(C10,X!2)),(C01,X(3)), 


: < C00+C1Q*XS0+C01*YSQ + 


1 (Cll ,X ( A) ) , (C20,X(5) ) , (C02,X<6) ) , 

2 (D00,X(7)),(D01,X(8)),{D02,X(9)) 

R2 = R *R 
R4=R2*R2 

: XSO=XI*XI 

X3= XSQ*X I 
X4=XSO*XSQ 
YSO= Y I *Y I 
Y3 = YSO*Y I 
Y4=YS0*YSQ 

P=GAMMA*XI+(OMGAM/R2>*X3+ 

1 XI* YSO* (R2-XSQ )* (DELTA+ ! R2-YSO) 

2 Cll *XS0*YSG+C20*X4+C02*Y4) ) + 

3 X3* ( R2-XSO) * (OMEGA+ ( R2-XSQ)* ( 000+ DO 1*X SO+D02*X 4 ) ) 
E=GAMMA*YI+(OMGAM/R2)*Y3+ 

1 YI*XSO* ( R2-YS0)* (DELTA+ ( R 2-XS 0 ) * ( C00+C10*YS0+C0 1*XSQ+ 

2 Cl 1*YS0*XSQ+C20*Y4+C02*X4) ) + 

3 V3* ( R 2- YSO ) * (OMEGA+ ( R2— Y S Q ) * (000+ DO 1* YSQ+002*Y A ) ) 

RETURN 
END 

CALCULATE THE AREA OF A TRIANGLE FORMED BY THREE CARTESIAN POINTS 
FUNCTION AREA (PI ,P2 ,P3) 

IMPLICIT REAL*8 ( A — H « 0 — Z ) 

DIMENSION PI (2 ) ,P2( 2 ) ,P3 (2 ) 

ASO=( PI ( 1) -P2 ( 1 1 )**2+IPl (2) -P2I2) )**2 
BSO= ( P2 ( 1 )-P3 ( 1 ) )*#2+ (P2(2) -P3 (2).)**2 
CS0*(P3(1)-P1(1))*#2+(P3(2)-P1(2)) **2 
B=DSQRT ( RSQ ) 

C=DSORT(CSQ) 

COS A=-( ASO-BSO-CSO) / (2.0D0*B*C) 

CP=B*COS A 
H=BS0-CP**2 

IF(H.LE,O.ODO) GO TO 100 
H= DS OR T (H) 

AREA=0, 500*C*H 
RETURN 

100 ARE A=O.ODO 
RF TURN 
FND 


00013700 

00013800 
00013900 
000 1 4000 
000 1 A 1 00 
00014200 
00014300 
00014400 
00014500 
00014600 
00014700 
00014800 
000 14900 
00015000 
00015100 
00015200 
00015300 
00015400 
00015500 
00015600 
00015700 
00015800 
00015900 
00016000 
00016100 
00016200 

00011900 
00012000 
00012100 
000 1 2200 
00012300 
00012400 
00012500 
00012600 
00012700 
00012800 
00012900 
00013000 
00013100 
00013200 
00013300 
00013400 
00013500 
00013600 
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E. 3 PROGRAM USED FOR EVALUATING COEFFICIENTS OF f*(£, 7?) 


MAIN 


TRjNSFnR»"T10N n R P .NcffoN H Iv C SlNI»IZTNG S ,HE dV,U>? OF 

A retransfgrmeo grid 


PROGRAMMER 


E. MICHAEL O'NEILL 
COMPUTER SCIENCES 
R7?R r.ni PRVTU.F R n 


CORPORATION 
. SILVER SPRING, 


MD. 


C 


C 

C 


C 

c 

c 

c. 

c 

c 

c 

c 


c 


c 


IMPLICIT REAL*R <A-H,0-Z) 

EXTERNAL FUNK 

COMMON / F ROOD/ ~*N FM AX » N F LA T »' JV AR Y , NE X TR A . 

COMMON /BILBD/^DIFF^VARj ^0^ 20? , N NT » MAT I MASK? 26 ) 

SFT DIMFNSIONAL CONSTANTS 

RS=1 .000 

R = 0. 5773502691896 DO 
RS0=0. 333333333333300 

R4=0. 11111 1 11111 1100 
INITIALIZE COEFFICIENTS 
DO 10 1=1,20 
V A R ( I 1=0.000 
VM I N ( T ) =-0.503 
VM A X ( I )= 0.5D3 

nv< I 1 = 0.100 

D E L M I N ( I 1=0.10-10 
10 MASK! 11=0 

NUMBER OF COEFFICIENTS 
N V = 1 3 

NUMBER OF ITERATIONS 
NFMAX=2000 
NT = 0 
M AT = 0 

X IM=R/ 5 .000 
Y IN = X IN 

STFPIT CALLS SUBROUTINE FUNK WHICH C A L CUL AT E S THE R M S 
DEVIATION OF THF RF TRANS FORMED GRID, STEPIT I J ^R AT ES U NT I L 
A SATISFACTORY RMS HAS BEEN REACHED OR UNTIL NFMAX ITERATIONS 

CALL STEP IT ( FUNK ) 


00000400 

00000450 

00000500 

00000550 

00000600 

00000700 

00000800 

00000900 

00001000 

00001100 


00001200 

00001300 

00001400 

00001500 

00001600 

00001700 

00001800 


00001900 

00002000 

00002100 

00002200 

00002300 

00002400 


00002500 


USE FINAL COEFFICIENTS TO CALCULATE RE TRANS FORMED GRID AND PRINT 



x=o.ono 

00002600 


00002700 


no 20 1 x= 1 , 6 

00002800 


y=o.ooo 
on 15 1 y= 1 , ix 

00002900 


DIRECT TRANSFORMATION 

00003000 


p=F0MAP (X ,Y 1 

E=FDMAP ( Y , X 1 

00003100 


INDIRECT TRANSFORMATION 

00003200 


XR=F IMAP ( P ,F 1 

00003300 


YR=F IMAP { F , P 1 

DIFF = ( X-XB 1 *4 2+ ( Y-YB) 

00003400 


00003500 


WRITF(6»12) 

00003600 


WRITE! 1,12) X,Y,P,E,XB,YB,DIFF 

FORMAT! I X , Y , P , F , XB , YB,D IFF ',7F8.51 

00003700 

12 

00003800 

15 

Y=Y+Y IN 

00003900 

20 

X=X+XIN 

00004000 

STOP 

END 

00004100 
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SUBROUTINE FUNK 


00007000 


C 

C THIS SUBROUTINE CALCULATES THF SUM OF THE SQUARES OF THE 

C DEVIATIONS OF THE RETRANSFORMED GRID 

C 

IMPLICIT REALMS <A-H,0-Z) 

COMMON/ R I LBO / 0IFF,VAR(20) 

COMMON RS*R»RS0,R4,XIN,YIN 
niFF=0.000 

x=o.ono 

00 10 IX=1,6 

y=o.odo 

00 5 I Y= 1 , IX 

C DIRECT TRANSFORMATION 

P=FDMAP (X, Y) 

E = FDMAP ( Y ,X ) 

C INVERSE TRANSFORMATION 

XR=F IMAP ( P,E ) 

YB=FIMAP(E,P) 

D I F F = D IFF + I X-XB 1**2+ ( Y-YB ) **2 
5 Y=Y+YIM 
10 X=X+XIN 
RETURN 
END 


00007100 

00007200 

00007100 

00007400 

00007500 

00007600 

00007700 

00007800 

00007900 

00008000 

00008100 

00008200 

00008B00 

00008400 

00008500 

00008600 

00008700 
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00005500 


FI) NC T I ON FDM&P ( A , B ) 

THIS FUNCT ION EVALUATES THE D I RFC T T R AN SF ORM AT I ON AT A,B USING 
A PREDETERMINED SET OF COEFFICIENTS 

IMPLICIT RF A L *8 (A-H.O-Z) 

8ATA n GAMMA^MGA^o:72360125A5882DO,O.2763987454A18DO/ 

DATA DELTA .OMEGA /0 .7 904 8 6449120800,-0. 1 2 2 544 1 4 R 7984D 1 / 

SaTA Doloi ^D2/0. 14833 1292 D1 ,0. 1 1 1997 260D1 » 0.60 51 5382101/ 

ASQ= A*A 

A3=A*AS0 

A4=AS0-:=AS0 

BS0=R4fi 

B4=RS0#BSQ 

RM A= R S 0- A SO 

iFD MAP-GAMMA *^gKGAM/RS0 ; f3+ + Di # As O+D2 *A4n + 

I f C00+ C 104 ASO+CO 1*BS0+C 1 14 ASQ4 BSQ+C204 A4+C024B4 ) ) 

RETURN 

END 

FUNCTION F I M AP ( A , R ) 

tut q FUNCTION EVALUATES THE INVERSE TRANSFORMATION USING A 
SET OF COEFFICIENTS INPUTTED THROUGH COMMON/ B I L B 0/ 


IMPLICIT REAL 
COMMON RS,R,R! 
COMMOM/R I LRO/ 
EQUIVALENCE If 
1 (Di 

1 <' 
1 1 1 
1 H 

DATA GAMMA, OMi 
ASO= A4A 
A3= A*A SO 
A4=AS0*ASQ 
R S 0= B * R 
B4=BS0*BS0 
RMA=RSQ-ASO 
RMA=RS0-AS0 
F IMAP=GAMMA*A 
1 A3*RMA>*( 

7 A *6 S04RM 

3 (C00+C10 

RETURN 
END 


48 ( A — H ,0-Z ) 

SQ.R4 ,X IN, Y IN 
D IF ,VAR ( 20) , 

OELTAtVARI 1) ) , (OMEGA, VAR ( 2 ) , 

SA 1 v i bI 1 * > ! : i & 2 i aS o 1 ! ! : < SSf lXS 1 i I? , • , 

IGAM/O. 13819765978860 1 , -0. 38 197659788 5500/ 


0MFGA + 0MEG1*(RS0-BS0)+RMA4 ( DO + D 1 A S 0+ D 2 m A 4 ) ) + 
IA4 (DELTA + OFI T1 4RM A+ ( RSO-BSO ) =■= 
i* AS0+C01*B SO+C 1 1 4 AS 04 B SQ+ C2 04 A4+ GO 24 B4 ) ) 


00005600 

00005700 

00005710 

00005720 

00005730 

00005740 

00005750 

00005800 

00005900 

00006000 

00006100 

00006200 

00006300 

00006400 

00006500 

00006600 

00006700 

00006800 

00006900 

00004200 


00004300 

00004400 

00004500 

00004525 

00004550 

00004600 

00004650 

00004700 

00004805 

00004810 

00004815 

00004820 

00004825 

00004830 

00004835 

00004900 

00005000 

00005100 

00005150 

00005200 

00005300 

00005400 
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E.4 MINIMIZER SUBROUTINE USED BY THE PROGRAMS OF SECTIONS E.2 
AND E.3 


SUBROUTINE STEP IT *FUNK< 

IMPLICIT REAL*8 (A-H.O-Z) 

360 STEPIT, MODIFIED BY J. BENSINGER, U. PENNSYLVANIA, 1971 

STEPIT FINDS LOCAL MINIMA OF A SMOOTH FUNCTION OF SEVERAL PARAMETERS, 

3STEPIT IS A PHLEGMATIC METHOD OF SOLVING A PROBLEM. 3 

— J. H. BURRILL, JR., 3360 STEPIT - A USER3S MANUAL3 

THIS SOURCE DECK AND A LONG WRITE-UP ARE AVAILABLE FROM .... 

QUANTUM CHEMISTRY PROGRAM EXCHANGE 
DEPT. OF CHEMISTRY, INDIANA UNIVERSITY 
BLOOMINGTON, INDIANA 47401 


4 4 4*4 
FUNK 


NV 

NTRACE 


MATRIX 


CHI SO 
MASK* J< 
X*J< 

XMAX* J< 
XMIN%J< 
0ELTAX*J< 
DEL M I N*J< 

ERR* J ,K< 

NFMAX 
NFL AT 

J VARY 


4444444444444444444444444444 

— THE NAME OF THE SUBROUTINE THAT COMPUTES CHISO 
GIVEN X*K,X*2C,...,X*NV< 

— THE NUMBER OF PARAMETERS, X 

— #0 FOR NORMAL OUTPUT, #tl FOR TRACE OUTPUT, 

#-l FOR NO OUTPUT 

— #0 FOR NO ERROR CALCULATION, #100£M FOR ERROR 
CALCULATION USING STEPS 1044M TIMES LARGER 
THAN THE LAST STEPS USED IN THE MINIMIZATION 

— THE VALUE OF THE FUNCTION TO BE MINIMIZED 
— NONZERO IF X*J< IS TO BE HELD FIXED 

— THE J-TH PARAMETER 

— THE UPPER LIMIT ON X*J< 

— THE LOWER LIMIT ON X*J< 

THE INITIAL STEP SIZE FOR X*J< 

THE LOWER LIMIT *CONVERGENCE TOLERANCE< ON THE 
STEP SIZE FOR X*J< 

RETURNS THE ERROR MATRIX *ALSO USED FOR SCRATCH 
STORAGES 

THE MAXIMUM NUMBER OF FUNCTION COMPUTATIONS 
NONZERO IF THE SEARCH IS TO TERMINATE WHEN ALL 
TRIAL STEPS GIVE IDENTICAL FUNCTION VALUES 
STEPIT SETS JV ARY NONZERO IF X*JVARY< IS THE ONLY 
X*J< THAT HAS CHANGED SINCE THE LAST CALL TO 
FUNK 

NEXTRA — USED BY SUBROUTINE SIMPLEX BUT NOT BY STEPIT 

*#44**44444444444444444444444444444 

THF DIMENSIONS OF ALL VECTORS AND MATRICES ARE NV, EXCEPT FOR .... 

E RR*NV ,MAX*NV .MOSQUE <<,SEC0ND*2 , 2<,X0SCSNV , MO SQUE< , C HI OSC*MOS DUES 

TO RFDUCE STORAGE TO A MINIMUM, SET M0S0UE#0, REDIMENSION ER R*1 , 1< » 
XOSCSl.K, SALVO* 1< , AND CHI0SC*1<, AND DELETE THE OSCILLATION 
SEARCH AND ERROR COMPUTATION SECTIONS *SEE COMMENT CARDS BELOWC. 

OR, USE SUBROUTINE STP, WHICH HAS THESE DELETIONS PLUS DELETION OF 
THE COLINEARITY CHECK. 

444444*4444444444444444444444444444 

/ 

COMMON /BILBO/ CH I SQ ,X < 2 0 > , XM IN ( 20 ! , X MAX ( 20 > , DEL TA X ( 20 ) , 

£ DELMIN (20 ) , E RR < 20 , 20 i , N V ,NTR ACE , MA TR I X , MA SK ( 20 ) 

COMMON /FRODO/ NFM AX ,NFL AT , J VAR Y , NEX TR A 

DIMENSION VEC*2O<,TRIAL*20<,XSAVE*2O<,CHI*2O<,DX*2O<,SECONDX2,2< 
DIMENS ION OL D VEC *20< ,SALVO*2 0<,XOSC?2O, 1 5< , CH I OSC* 1 5< , J FL AT*20< 

TO RUN ON THE CDC 3600, THE FOLLOWING DATA STATEMENT MUST BE 
ACTIVATED AND SUBROUTINE SSWTCH MUST BE SUPPLIED. 

TO RUN ON OTHER COMPUTERS, A BLOCK DATA SUBPROGRAM MUST BE SUPPLIED 
TO SET NFMAX AND NFLAT. 

DATA *NFMAX#1000000<,*NFLAT#1< 

4*444444444444444444444444444444444 

SET FIXED QUANTITIES .... 

KW ... STANDARD OUTPUT UNIT NUMBER 
IBM 709/7090/7 09 A .... 

KW #6 

CDC 3600/6400 .... 

K W #6 1 

KTYPE ...- CONSOLE TYPEWRITER UNIT • 


GRH03140 

GRH03150 

GRH03160 

GRH03180 

GRH031A0 
GRH031B0 
GRH031CO 
GRH031D0 
GRHO31E0 
GRH031F0 
GRH03200 
GRH032 1 0 
GRH03220 
GRH03230 
GRH03240 
GRH032 50 
GRH03260 
GRH032 70 
GRH03280 
GRH03290 
GRH032A0 
GRH032B0 
GRH032C0 
GRH032D0 
GRH032E0 
GRH032F0 
GRH03300 
GRH033 1 0 
GRH03320 
GRH03330 
GRH03340 
GRH033 50 
GRH03360 
GRH03370 
GRH03380 
GRH03390 
GRH033A0 
GRH033B0 
GRH033C0 
GRH033D0 
GRHO33E0 
GRHO33F0 
GRH03400 
GRH03410 
GRH03420 
GRH03430 
GRH03440 
GRH03450 
GRH03460 
GRH0347 0 
GRH03480 
GRH03490 
GRH034A0 
GRH034B0 
GRH034C0 
GRHO34D0 
GRH034E0 
GRH034F0 
GRH03500 
GRH03510 
GRH03520 
GRH03530 
GRH03540 
GRHO3550 
GRH03560 
GRH03570 
GRH03580 
GRH03590 
GRH035A0 
GRH035B0 
GRH035C0 
GRH035D0 
GRH035E0 
GRH035F0 
GRH03600 
GRH0361 0 
GRH03620 
GRH03630 
GRH03640 
GRH03650 
GRH03660 
GRH03670 
GRH03680 
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K TYPE #6 


SIGNIFW2 ,E8 

OOUBLE PRECISION 
S IGNI F#1 .016 

SIGNIF#1.E10 

SIGNIF#1.E16 


HUGE # 1 . E 37 
HUGE # 1 * D75 
HUGE # 1 . E307 

NVMAX#20 

MOSOUE#1 5 
RATI0#10.0 
OIL IN#0.99 

Nf,nMP*5 

ACK#2.0 

F ACUP #4 . 
NSTPMX#5 
NSSW#6 

*$**#*** 


COC 3600 .... 

SIGNIF ... LARGEST NUMBER SUCH THAT 
XSIGNIF£1.<-SIGNIF IS NONZERO 

IBM 709/7090/7096 .... 

I BM 360 

COC 3600 .... 

COC 6600 .... 

HUGE ... THE LARGEST FLOATING POINT 
NUMBER 

IBM 709/7090/7096 .... 

IBM 360 

COC 3600/6600 .... 


NVMAX ., 

MOSQUE ... 

SEARCH 


MAXIMUM VALUE OF NV 
MAXIMUM DEPTH OF OSCILLATION 

RATIO ... RATIO OF STEP SIZE DECREASE 

COLIN ... COLINEARITY TOLERANCE 

NCOHP ... MAXIMUM NUMBER OF CYCLES 
WITHOUT ATTEMPTING A GIANT STFP 

ACK ... RATIO OF STFP SIZE INCRFASF 

FACUP ... IF MORE THAN FA CUP STEPS ARE 
TAKEN, THE STEP SIZE IS I NCREASFD 

NSTPMX ... L0G2XMAXIMUM NUMBER OF STEPS< 

NS SW ... TERMINATION SENSE SWITCH NUMBER 


CHFCK SOMF INPUT QUANTITIES, AND SET THEM TO STANDARD VALUES IF 
DESIRFD. 

DECOMMENT FOR SSW OPTION .... 

CALL SSWTCH XNSSW, JUMPS 
IFXJUMP-1S10 ,10 ,60 

ONLY USAGE OF THE CONSOLE TYPEWRITER.... 

10 WRITE XKTYPE,20< NSSW 

20 F0RMATX23H TURN OFF SENSE SWITCH I 2< 

30 CALL SSWTCH XNSSW, JUMPS 
IFXJUMP-1<30,30,60 

60 IFXNV<630,630,50 
50 IFXNCOMPS60 ,60 ,70 
60 NCOMP U 1 
70 J VARY #0 


N ACT IV . 

NAC T I V #0 
DO 230 I # 1 , NV 
I FXMASKXI <<2 30,80,2 30 
80 DFLX#DELTAXXI< 

I FXDELX*XXI <<90,100,1 00 
90 DE L X#— DE L X 
100 XPLUS#XXI<£DELX 

IFXXPLUS-XXKS160, 110,160 
110 IF? XX I«1 30, 120, 130 
120 DELTAXXI<#0.01 
GO TO 160 

130 DELTAXXI<#0.01*XXI < 

160 rFXDELMINXI«160,150 , 170 
150 OELMINXI<#DELTAXXI </SIGNIF 
I FXD EL MIN XI <<16 0,170 ,170 
160 DELMINXI<#-DELMINSI< 

170 I FXXMAXXK-XMINXI <<180,180, 19 0 
180 XMAXXI <#HUGE 


NUMBER OF ACTIVE XXJ< 


GRH03690 

GRH036A0 

GRH036B0 

GRH036C0 

GRH036D0 

GRH036E0 

GRHO36F0 

GRH03700 

GRH0371 0 

GRH03720 

GRH03730 

GRH03760 

GRH03750 

GRH03760 

GRH03770 

GRH03780 

GRH03790 

GRH037A0 

GRH037B0 

GRH037C0 

GRH037D0 

GRH037F0 

GRH037F0 

GRH03800 

GRH03810 

GRH03820 

GRH03830 

GRH03860 

GRH03850 

GRH03860 

GRH03870 

GRH03880 

GRH03890 

GRH038A0 

GRH038B0 

GRHO38C0 

GRHO38D0 

GRH038E0 

GRHO38F0 

GRH03900 

GRH0391 0 

GRH03920 

GRH03930 

GRH03960 

GRH03950 

GRH03960 

GRH03970 

GRH03980 

GRH03990 

GRH039A0 

GRH039B0 

GRH039C0 

GRHG39D0 

GRH039E0 

GRH039F0 

GRH03A00 

GRH03A 10 

GRH03A20 

GRH03A30 

GRH03A60 

GRH03A50 

GRH03A60 

GRH03A70 

GRH03A80 

GRH03A90 

GRH03AA0 

GRH03AB0 

GRH03B60 

GRH03B50 

GRH03B60 

GRH03B70 

GRH03B80 

GRH03B90 

GRH03BA0 

GRH03BB0 

GRH03BC0 

GRH03BD0 

GRH03BEO 

GRHO3BF0 

GRHO3C0O 

GRH03C 1 0 

GRH03C20 

GRHO3C30 

GRH03C60 
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XM IN? I <#— HUGE 
190 NACTIV#NACTIV£1 

X?I <#AMAX1?XMIN?K,AMIN1?XMAX?I<,X?I<« 
IF?X?I<-XMAX?I«210,210,200 
200 X? I <#XMAX ?I < 

210 I F?X ? I <— XM I N?I <<220,230 *230 
220 X? I <#XM IN?I < 

230 CONTINUE 

COMPUTE C0MPAR. THE PROBABILITY THAT THE COSINE OF THE ANGLE 
BETWEEN TWO RANDOM DIRECTIONS EXCEEDS COMPAR IS APPROXIMATELY 
? 1-COL I N</ 2 . 

COMP AR#0 • 0 

IF?NACTIV— 1 <240 ,300,260 
240 DO 250 J#1 ,NV 
250 MA SK? J <#0 
GO TO 50 
260 A# NACT I V 

SUB#2.0/?A-1.0< 

P#2 .0**1 .0/SQRT ?A</?1 .0-0.5 **SUB<-1.0< 

P#2.0*?l .0/DSQRT?A</?l . 0-0.5** SUB<-1 .0< 

C0MPAR#?1 .-?1 .-C0LIN<**SUB<*?1. £P*?1.-C0LI N<< 

C0MPAR#MIN1F?.999,ABSF?C0MPAR<< 

IF?C0MPAR<270 ,280 ,280 
270 C0MPAR#-C0MPAR 
280 IF?CDMPAR-.999<300,300,290 
290 COMPAR#. 999 

300 IF?NTRACE<390,310,310 


310 WRITE? r\w,^4us uRnu^ir^u 

320 FORM A T?92H1 ENTER SUBROUTINE STEPIT. ORIGINAL VERSION BY CHANDL ERGRH03E40 


KW , 32 0< 


GRH03C50 
GRH03C60 
GRH03C70 
GRH03C80 
GRH03C90 
GRH03CA0 
GRH03CB0 
GRH03CC0 
GRH03CD0 
GRH03CE0 
GRH03CF0 
GRH03D00 
GRH03D 1 0 
GRH03D20 
GRH03D30 
GRH03D40 
GRH03D50 
GRH03D60 
GRH03D70 
GRH03D80 
GRH03D90 
GRH03DA0 
GRH03DB0 
GRH03DC0 
GRHO3DD0 
GRH03DE0 
GRH03DF0 
GRH03E00 
GRH03E10 
GRH03E20 
GRH03E30 


*, PHYSICS DEPT., INDIANA UN I VF R S IT Y . // 19H INITIAL VALUES. .../< 


WRITE? 

330 F0RMAT?/10H MASK 
up |TFS 

340 FORMAT?/ 10H X 
WR I TF ? 

350 FORMAT?/ 10H XMAX 
WR I TE ? 

360 FORMAT ?/ 1 OH XMIN 
WRITE? 


KW ,330<?MASK?J< , J#1 ,NV< 

# 9?I6 ,6X</?4X9I 12« 

KW ,340<?X?J< , J#1 , NV < 

# 9E12 ,4/?10X 9E12.4« 

KW ,350<?XMAX?J< , J#1 ,NV< 

# 9E12.4/?10X 9 E 1 2 . 4<< 

KW ,360<?XM I N?J< , J#1 ,NV< 

# 9E12.4/?10X 9 E 1 2 • 4<< 

KW ,37 0<?DELTAX?J<, J#1 ,NV< 

370 FORMAT ?/ 10H DELTAX # 9F12.4/?10X 9E12.4<< 

WRI TF? KW,380<?DELMIN?J<, J#1 ,NV< 

380 FORMAT ?/ 10H DELMIN # 9E12.4/?10X 9E12.4<< 

390 CALL FUNK 

CHISAV#CHISO 
CALL FUNK 

NF ... NUMBER OF FUNCTION CALLS 

NF #2 

IF?CHISO-CH ISAV<400 ,420 ,400 
400 WRITE ?KW , 410< CH I SAV ,CH I SO ,NF 
410 F0RMAT?///3?/59H WARNING.... CH I SO IS NOT A 
* OF X?JC. < A5X 9HCHISAV # E22.14.5X 8HCH I SO 


GRH03E50 
GRH03E60 
GRHD3E70 
GRH03E80 
GRH03E90 
GRHD3EA0 
GRH03EB0 
GRH03EC0 
GRH03 EDO 
GRH03EE0 
GRH03EF0 
GRH03F00 
GRH03F10 
GRH03F20 
GRH03F30 
GRH03F40 
GRH03F50 
GRH03F60 
GRH03F70 


GRH03F80 

REPRODUCIBLE FUNCT 1 0NGRH03 F90 
# E22.14,5X 5HN F # I5<GRHO3FA0 
GRH03FB0 

JOCK ... SWITCH USED IN SETTING JVARY GRH03FC0 

420 J0CK#1 GRH03 EDO 

IF?NTRACE<450»430,430 GRH03FE0 

430 WRITE? KW , 440<N V ,NACTIV, MATRIX , NCO MP,NFMAX,NFLAT, GRHO3FF0 

* RATIO, ACK, COLIN, COMPAR, CHISO GRH04000 

440 FORMAT ?/ / 1H I3.11H V A R I A B L E S , 1 3 , 8H ACT I V E . 1 OX 8HM ATR I X #14, 10 X7HNC0GRH040 1 0 

*MP #12, 10X7HNFMAX # I 8 , 1 0X7HNFL A T #12// GRH04020 

* 8H RATIO #F5 . 1 , 1 OX 5H ACK #F 5 . 1 , 10X 7H CO L I N #F6 . 3 , 10X 8HC OM PA RGRH04030 

* #F6.3///8H CHISO #E15.8///23H BEGIN MINIMIZATION.... < GRH04040 

450 IF?NV<2120, 2120, 460 GRH04050 

460 IF?NTRACE<490,490,470 GRH04060 

470 WRITE? KW ,480< GRH04070 

480 FORMAT?//60?1X1H*<//10X37HTRACE MAP OF THE MINIMIZATION P ROC E S S/ / <GRH04080 

GRH04090 
GRH040A0 

VEC? J< ... CURRENT VECTOR OF NUMBER OF GRH040B0 
STEPS IN X?J< GRH040C0 

GRH040DO 

DX?J< ... CURRENT STEP SIZE FOR X?J< GRH040EO 

GRH040F0 
GRHO410O 
GRH041 10 

NOSC ... CURRENT DEPTH OF THE OSCILLATION GRH04120 
INFORMATION GRH04130 

GRH04 1 40 

KW I T ... TERMINATION SWITCH GRH04150 

GRH04160 
GRH04170 

********************* GRH04180 


490 DO 500 I # 1 , NV 


VEC?I <#0. 

500 DX?I<#DELTAX?I< 
CHI OLD #CH I SQ 


N0SC#0 
KW I T #0 


CHIOLD ... BEST PREVIOUS VALUE OF CHISO 


***####*#***** 
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VARY THE PARAMETERS ONE AT A TIME. 

THIS IS THE STARTING POINT USED EACH TIME THE STEP SIZE IS REDUCED 
OR A SUCCESSFUL GIANT STEP IS COMPLETED. 


510 NCIRCOO 


N i I P HO 


NC IRC ... NUMBER OF CONSECUTIVE XSJ< 
WITHOUT SIZEABLE CHANGES 


NZIP 


... NUMBER OF CONSECUTIVE CYCLES 
WITHOUT A GIANT STEP 


MAIN DO LOOP FOR CYCLING THROUGH THE VARIABLES. 
FIRST TRIAL STEP WITH EACH VARIABLE IS SEPARATE. 


HACK 


520 NACKHO 
DO 1750 


... NUMBER OF ACTIVE 
THROUGH 


X*J< CYCLED 


I i'll »NV 
FLAT SI <s¥0 


J F L A T % J < ... NONZERO IF CHANGING XSJ< DID 
NOT CHANGE CHI SO 


OLDVECSJ< 


DLD VECTOR OF NUMBER OF 


OLDVECSI <>1VECSI < 

VEC2I <#0.0 

TR I AL SI <#0 . 0 
I FSMASKSI <<530,540,530 
530 VECSI<#-0.0 
JFLATSI <#1 
GO TO 1730 
540 NACICSNACKSl 
S A V6 #X SI < 


STEPS IN XSJ< 

TRIALSJ< ... CHANGE IN XSJ< 


CHECK THAI DXSI< IS NOT NEGLIGIBLE. 


DE LX#DX % I < 

IFSDELX*SAVE<550»560 ,560 
550 DELX#-DELX 
560 XPLUS#SAVE£DELX 

IFSXPLUS-SAVE<580, 570,580 
570 JFLATSI <#2 
GO TO 770 

ADX#ABSSDXSI « 

580 ADX #DXSI < 

IFSADX<590, 600,600 
590 ADX # — ADX 

CTCD V.t/ 

600 XSI <#SAVEEDXSI < 

JVARY «0 

IFSJGCK<620,620»610 
610 JOCK #0 
JVARY #1 
620 NFL AG #1 

I ESXST<-XM I NS I <<640, 6 30 ,6 30 
6 30 IFSXS I <-XMAXSI <<6 50 , 6 50 ,6 4 0 
640 NFl AG#MFLAG£3 
GO TO 670 
650 CALL FUNK 
N F #N F £ 1 
JVARY #1 

SAVE OLD VALUE OF CH I SO FOR INTERPOLATION. 

CH I ME #CH I SQ 

I FSCH I SQ—CH I0LD<810 ,660 ,670 
660 NFL AG#NFLAG£ 1 

STEP HI< THE OTHER WAY. 

670 XSI <#SAVE-DXSI< 

IFSXS I <-XMINSI <<7 80 ,6 80 ,6 8 0 
680 IFSXSK“XMAXSI«690,690,780 
690 CALL FUNK 
N F #N F £ 1 
JVARYKI 

I FSECH I SQ— CHI0LD<800 ,700 ,7 10 
700 NFL AG KNFL AGE 1 
710 IFSNFLAG— 3<720,760,780 

PERFORM PARABOLIC INTERPOLATION. 

CHECK FOR ZERO DENOMINATOR, ETC. 

720 IFS?35CHISQ~CHIME<>i'SCHIME-2.*CHIOLD£CHlSQ«7 30,780,7 30 
730 TRIALSI<#.5*DXSI <*SCH ISQ-CH IME</SCH IM E-2 . *CH I OL D£ CHI SQ< 

VECSI<#TRI ALSK/ADX 
X%I<#SAVE£TRIAL*I< 

CALL FUNK 
NFgNFE 1 

I FSCHISQ--CHIOLD<740,750 ,750 


GRH04190 
GRH041A0 
GRH041B0 
GRH041C0 
GRH04100 
GRHO41E0 
GRH041F0 
GRH04200 
GRH04210 
GRH04220 
GRH04230 
GRH04240 
GRH04250 
GRH04260 
GRH04270 
GRH04280 
GRH04290 
GRH042A0 
GRH042B0 
GRH042C0 
GRH042D0 
GRH042E0 
GRH042F0 
GRH04300 
GRH04310 
GRH04320 
GRH04330 
GRH04340 
GRH043 50 
GRH04360 
GRH04370 
GRH04380 
GRH04390 
GRH043A0 
GRHO43B0 
GRH043C0 
GRH043D0 
GRH043F0 
GRHO43F0 
GRH04400 
GRH04410 
GRH04420 
GRH04430 
GRH04440 
GRH044 50 
GRHD4460 
GRH04470 
GRH04480 
GRH04490 
GRH044A0 
GRH044B0 
GRH044C0 
GRH044D0 
GRH044E0 
GRH044F0 
GRH04500 
GRH0451 0 
GRH04520 
GRH04530 
GRH04540 
GRH04550 
GRH04560 
GRH04570 
GRH04580 
GRH04590 
GRH045A0 
GRHO45B0 
GRH045C0 
GRH04 5D0 
GRH045E0 
GRHO45F0 
GRH04600 
GRH04610 
GRH04620 
GRH04630 
GRH04640 
GRH04650 
GRH04660 
GRH04670 
GRH04680 
GRH04690 
GRH046A0 
GRH046B0 
GRH046C0 
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740 CH IOLD #CH ISO 
JOCK# 1 
GO TO 790 
750 TRIAL SI <#0.0 
VEC SI <#0 .0 
GO TO 780 
760 JFLATSI <#1 
770 VECSI <#-0. 

780 XS I <#SAVE 
790 NC I RC #NC I RC£ 1 

IFSNCIRC-NACTIV<910,1820, 1820 

C RESET OXS I < FOR MORE EFFICIENT OPERATION. 

800 DXSI<#-DXSI< 

A LOWER VALUE OF CHISQ HAS BEEN FOUND. STEP, DOUBLE THE STEP SIZE, 
AND REPEAT AS LONG AS CHISQ CONTINUES TO DECREASE. 

810 NC IRC#0 
DEL#DXSI< 

NSTP#0 

820 CH IME#CHIOLD 
CH I OLD#CH I SQ 

VECSI<#VECSI<£DEL/ABS SDXSI« 

VECSI<#VECSI<6DEL/DABSSDXSI« 

TRIAL SI <#TRIALSK£DEL 
NSTP#NS TP £ 1 

IFSNSTP— NSTPMX<830,890,890 
830 OEL #ACK*DEL 
SAVE #X SI < 

XSI <#SAVE£DEL 

IFSXSI<-XMINSI<<900 ,840 ,840 
840 I FSXSI <-XMAXSI <<850 ,850 ,900 
850 CALL FUNK 
N F #N F £ 1 

IFSCHISO-CHIOLO<820,86O ,860 

PERFORM PARABOLIC INTERPOLATION. 

860 DEN0M#ACK*CHIME-SACK£1.<*CHI0LD£CHI SQ 
I FSDEN0M<87 0,900 ,870 

870 C INDER#S.5/ACK<*SACK**2*CHIME-SACK**2-1.<*CHI0LD-CHI SO</DENOM 
XSI <#SAVE£CINDER*OEL 
CALL FUNK 
N F #N F £ 1 

I FSCH I SQ-CH IQLD<880,900 ,900 
880 CH I OL D#CH I SQ 

TRIALSI<#TRIALSI<SCINDER*DEL 
VECSI <# VECSI <£C I NDER*DEL/ADX 
890 JOCK #1 

GO TO 910 
900 XS I <#S A VE 

DO NOT INCREASE THE STEP SIZE PREMATURELY. 

910 IFSNZIP<920 ,920 .940 
920 IFSNCOMP-K930, 930, 1720 
930 IFSNACK-K1720, 1720, 940 

» AVEC#AB SS VECSI << 

940 AVEC#VECSI< 

IFSAVEC<950, 96 0,960 
980 AVFC*-AVFC 

960 TFSAVFC— FAC UP< 1030, 970, 970 
970 DXSI<#ACK*ADX 

VECSI <#VECSI</ACK 
OLDVECSI <#0L 0 VECS I </ ACK 
I FSNOSC< 1000, 1000,980 
980 DO 990 J # 1 ,NOSC 
990 ERRSI ,J<#ERRSI ,J</ACK 
1000 IFSNTRACE<1030, 1030 , 1010 
1010 WRITES KW,1020<I , DX S I < 

1020 F0RMATS10H STEP SIZEI3.14H INCREASED TO E12.5< 

*«❖##*****##*#**#####*#«********#«# 

IF POSSIBLE STEP ALONG A RESULTANT DIRECTION. 

CHECK THE COLINEARITY OF VEC AND OLDVEC. 

1030 SUM0#0 . 0 
SUMV#0 .0 
DO 1040 J# 1 ,NV 
SUM0#SUM0£0LDVECSJ<**2 
1040 SUM V# SUM V£VECSJ<*42 

I FSSUMO*SUMV <1720 ,1720,1050 
C 10 50 SUMO#SQRT SSUMO< 

1050 SUMO#DSQRTSSUMO< 

C SUM V#SQR T SSUMV< 

SUMV#DSQRTSSUMV< 

COS I NE#0 . 0 


GRHO46D0 
GRH046E0 
GRH046F0 
GRH04700 
GRH04710 
GRH04720 
GRH04730 
GRH04740 
GRH04750 
GRH04760 
GRH04770 
GRH04780 
GRH04790 
GRH047A0 
GRH047B0 
GRH047C0 
GRH047D0 
GRH047E0 
GRH047F0 
GRH04800 
GRH04810 
GRH04820 
GRH04830 
GRH04840 
GRH04850 
GRH04860 
GRH04870 
GRH04880 
GRH04890 
GRH048A0 
GRH048B0 
GRH048C0 
GRH048D0 
GRH048E0 
GRH048F0 
GRH04900 
GRH04910 
GRH04920 
GRH04930 
GRH04940 
GRH04950 
GRH04960 
GRH04970 
GRH04980 
GRH04990 
GRH049A0 
GRH049B0 
GRH049C0 
GRH049D0 
GRH049E0 
GRH049F0 
GRH04A00 
GRH04A 10 
GRH04A20 
GRH04A30 
GRH04A40 
GRH04A50 
GRH04A60 
GRH04A70 
GRH04A80 
GRH04A90 
GRH04AA0 
GRH04AB0 
GRH04AC0 
GRH04AD0 
GRH04AE0 
GRH04AF0 
GRH04B00 
GRH04B 1 0 
GRH04B20 
GRH04B30 
GRH04B40 
GRH04B50 
GRH04B60 
GRH04B70 
GRH04B80 
GRH04B9D 
GRH04BA0 
GRH04BB0 
GRH04BC0 
GRH04BD0 
GRH04BE0 
GRH04BF0 
GRH04C00 
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DO 1060 J#1 ,NV 

1060 C0SINE#C0SINE£S0LDVECSJ</SUM0<*SVECSJ</SUMV< 
IFSNZIP-NC0MP<1070.1100 , 1 100 
1070 IFSNZIP-K1720, 1080, 1090 
1080 I F SNAC K-NAC TIV<1720,1090, 1090 
1090 IFSCOSINE-COMPAR<1720, 1 100, 1100 
1100 IFSVEC*I<<1110, 1720, 1110 
1110 NON ZER #0 

00 1130 J #1 »NV 
I F*VEC«J<< 1120, 1130,1120 
1120 N0NZER#N0NZER£1 
1130 CONTINUE 

IFSNONZER-2<1720,1140,1140 
SIMON SAYS, TAKE AS MANY GIANT STEPS AS POSSIBLE... 

NG I ANT ... NUMBER OF GIANT STEPS 


1140 NG I ANT#0 

NTR Y#0 

NR E TRY #0 
KL#1 


NT R Y ... NONZERO IF A GIGANTIC 

^OSCILLATIONS STEP IS BEING TRIED 

NRETRY ... NUMBER OF OSCILLATION PERIODS 
YET TO BE TESTED 

KL ... POINTER FOR OSCILLATION CHECK 
STORE OSCILLATION INFORMATION. 


N0SC#N0SC£ 1 
IFSNOS C-M OS QU £<1180,1180, 1150 
1150 N0SC#M0SQUE 

I FSNOSC- K1300.ll 80, 1160 
1160 DO 1170 K#2 ,NOSC 

CHI0SCSK-1<#CHI0SC%K< 

DO 1170 J# 1 , NV 
XOSCSJ ,K-l<#XOSCSJ,K< 

1170 ERRSJ,K-1<#ERRSJ,K< 

1180 CONTINUE 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

TO DELETE THE OSCILLATION SEARCH SECTION, REMOVE ALL STATEMENTS 
BETWEEN HERE AND THE NEXT COMMENT CARD OF X 3 S. 

DO 1190 J#1 ,NV 
XOSCSJ,NOSC<#XSJ< 

1190 E RR SJ ,NOSC <#VEC SJ </ SUMV 
CHIOSCSNOSC<#CHIOLD 
I FSNOSC-3 <1300, 1200, 1200 

SEARCH FOR A PREVIOUS SUCCESSFUL GIANT STEP IN A DIRECTION MORF 
NEARLY PARALLEL TO THE DIRECTION OF THE PROPOSED STEP THAN WAS THE 
IMMEDIATELY PREVIOUS ONE. THIS MAY MEAN THAT THE DIRECTIONS OF THE 
GIANT STEPS OSCILLATE PERIODICALLY. TRY GIGANTIC SOSC I L L AT 10N< 

STEPS OF DECREASING PERIOD, THEN ORDINARY GIANT STEPS. 

1200 COXCOM #0.0 

DO 1210 J#1,NV 

1210 COXCOM #COXCUM£E RR SJ ,NOSC<*E RR * J ,N0 S C- 1< 

NAH #NOSC— 2 
1220 NTRY#0 

DO 1240 K#KL ,NAH 
NRETRY #NAH-K 
COS I NE#0 . 0 
DO 1230 J #1 ,NV 

1230 COSINE#COSINE£ERR?J,NOSC<*ERR*J,K< 

IF %COSINE-COXCOM<1240, 1240, 1250 
1240 CONTINUE 

GO TO 1300 
1250 NTRY# 1 
KL#K£1 

I F%NTR AC E< 1280, 128 0,1260 
1260 NT #NOS C— K 

WRITES KW,1270<NT, COXCOM, COSINE 

1270 F0RMATS/1X8H********5X26HG IGANT IC STEP WITH PERIOD 12, 

4 37H BEING ATTEMPTED. COXCOM, COSINE # 2E13.4< 

1280 DO 1290 J#1 ,NV 

SALVOSJ<#TRIALSJ< 

1290 TRIALSJ<#SXSJ<-XOSCSJ ,K«/ ACK 

CHIBAK#CHIOLD£XCHIOSCSK<-CHIOLD</ACK 
GO TO 1310 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 


GRH04C 1 0 
GRH04C20 
GRH04C30 
GRH04C40 
GRH04C 50 
GRH04C60 
GRH04C70 
GRH04C80 
GRH04C90 
GRH04CA0 
GRH04CB0 
GRH04CC0 
GRH04CD0 
GRH04CE0 
GRH04CF0 
GRHD4D00 
GRH04D10 
GRH04D20 
GRH04D30 
GRH04D40 
GRH04D50 
GRH04D60 
GRH04D70 
GRH04D80 
GRH04D90 
GRH04DA0 
GRH04DB0 
GRH04DCO 
GRH04DD0 
GRH04DE0 
GRH04DF0 
GRH04E00 
GRH04E10 
GRH04E20 
GRH04E30 
GRHD4E40 
GRH04E50 
GRH04E60 
GRH04E70 
GRH04E80 
GRHD4E90 
GRH04EA0 
GRH04EB0 
GRH04EC0 
GRH04ED0 
GRHD4EE0 
GRH04FF0 
GRH04F00 
GRH04F10 
GRH04F20 
GRH04F30 
GRH04F40 
GRH04F50 
GRH04F60 
GRH04F70 
GRH04F80 
GRH04F90 
GRH04FA0 
GRH04FB0 
GRH04FC0 
GRH04FD0 
GRH04FE0 
GRH04FF0 
GRH05000 
GRH05010 
GRH05020 
GRH05030 
GRH05040 
GRH05050 
GRH05060 
GRH05070 
GRH05080 
GRH05090 
GRH050A0 
GRHO50BO 
GRH050C0 
GRH05000 
GRHO50EO 
GRH050F0 
GRH05100 
GRH051 10 
GRH05120 
GRH05130 
GRH05140 
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PERFORM GIGANTIC OR GIANT STEPS. 


1300 CH1BAK#CHI?I< 

1310 00 1360 J#1,NV 
XSAVE?J<#X?J< 

TRI AL?J<#ACK*TRI AL?J< 

IF?MASK?J«136 0 ,1320 ,1360 
1320 X?J<#X?J<£TRI AL?J< 

C X? J<#AM AX ISAM IN 1?X?J<,XMAX?J«, XM IN?J« 

IF?X?J<-XMAX?J<<13A0,13A0, 1330 
1330 X?J<#XMAX?J< 

13 AO IF?X?J<-XMIN?J«13 50, 1360, 1360 
1350 X?J<#XMIN?J< 

1360 CONTINUE 
J0CK#0 
J VARY #0 
CALL FUNK 
NF #NF £ 1 

I F%CH IS O-CHI OLD <1370 , 1AA0.1AA0 
1370 CHI BAK#CH IOLD 
CH IOL D #CH I SO 
NGI ANT#NGI ANT£1 
IFXNTRACE <1310.1310,1380 
1380 I F?NGIANT— 1<1390, 1390, 1A20 
• l^on hrttf? «<•' , lAnn<r.Hi bak.nf , ?vec?j<, j#i, i < 

1A00 F0RMAT?//8H CHISO #F16.fl .8XAHNF * I 111 5X 16HN0 . OF STEPS # 9E11.3/ 

A ?2 1 X9E 1 1 . 3<< 

WRITES! KW,1A10<?XSAVE?J<, J#1 ,NV< 

1 A10 FORMA T S9H X SI <. . . ./ ? IX 9E 1 3 . 5 « 

1A20 WRITE? KW,1A30<CHISQ.NF,?XSJ<, J#1,NV< 

1 A30 F0RMATS/8H CHISO #E 16 .8 ,8X AHNF #I7/9H X SI <. . . . /? IX 9 E 1 3. 5<< 

GO TO 1310 

C DO NOT INTERPOLATE AFTER AN UNSUCCESSFUL 

C GIANT STEP. 

1 AAO IFSNGIANT <1530 , 1530 ,1A50 

C PERFORM PARABOLIC INTERPOLATION. 

1 A50 DEN0M#ACK*CHIBAK-?ACK£1.<*CHI0LD£CH I SQ 
IF?DEN0M<1A60,1520, 1A60 

1 A60 C INDER#S.5/ACK<*?ACK**2*CHIBAK-SACK**?~1.<*CHI0LD-CHI SO</DENOM 
DO 1510 J# 1 , NV 
IF?MASK?J<<1510 , 1A70 ,1510 
1A70 X?J<#XSAVESJ<£CINDER*TRI AL?J< 

C XSJ<#AMAX 1?AMIN1?X?J<,XMAXSJ«,XMINSJ<< 

I FSXSJ<-XMAXSJ<<1A90,1A90, 1A80 
1 A80 XSJ<#XMAX?J< 

1A90 I F?X?J<-XM I N? J<< 1500, 1510, 15 10 
1500 X? J <#XM I N ?J < 

1510 CONTINUE 
JOCK #0 
JVARY #0 
CALL FUNK 
NF#NF£1 

I F?CH I SQ-CH IOLO <1650 ,1520,1520 
1520 IF?NGIANT<1560, 1530, 1560 
1530 IF?NTRY<15A0,1560,15A0 
15A0 DO 1550 J#1,NV 

TRIAL?J<#SALV0?J< 

1550 X?J<#XSAVESJ< 

GO TO 1580 
1560 DO 1570 J#1,NV 

TRIAL ?J<#TRIALSJ</ACK 
1570 XSJ<#XSAVE?J< 

1580 IF?NTRACE<1610 , 1610 , 1590 

1590 WRITE? KW,1600<CHI0LD,NGIANT 

1600 FORMAT %/ 8H CHISQ #E15.8,7H AFTERI3 , 13H GIANT STEPS. < 

WRITE? KW, 1A10<?XSJ<, J#1 ,NV< 

1610 IF?NG IANT<1620 , 1620 , 1680 
1620 I F?NRETRY<1630, 1630, 1220 
1630 I F?NTR Y <16A0 , 1700 ,16A0 

ALL GIGANTIC STEPS WERE UNSUCCESSFUL. 

TRY A GIANT STEP. 

16 AO NTRY#0 

GO TO 1300 
C 

1650 CH IOL D #CH ISO 
JOCK# 1 

IF?NTRACE<1680 ,1680,1660 
1660 STEPS#NG I ANT 

STEPS#STEPS£CI NDER 

WRITE? KW,1670<CHIOLD, STEPS 

1670 FORMAT ?/8H CHISO #E15.8,7H AFT ERF6 . 1 , 1 3H GIANT STEPS. < 

WRITE? KW , 1A10<?X?J<, J#1 ,NV< 

1680 IF?NTRY<1690, 510, 1690 
1690 CONTINUE 


GRH051 50 
GRH05160 
GRH05170 
GRH05180 
GRH05190 
GRH051A0 
GRH051B0 
GRH051C0 
GRH051D0 
GRH051E0 
GRH051F0 
GRH05200 
GRH05210 
GRH05220 
GRH052 30 
GRH052A0 
GRH05250 
GRH05260 
GRH05270 
GRH05280 
GRH05290 
GRH052A0 
GRH052B0 
GRH052C0 
GRH052D0 
GRH052F0 
GRH052F0 
GRH05300 
GRH05310 
GRH05320 
GRH05330 
GRH053A0 
GRH05350 
GRH05360 
GRH05370 
GRH05380 
GRH05390 
GRH053A0 
GRH053B0 
GRHO53C0 
GRH053D0 
GRHO53F0 
GRH053F0 
GRHO5A00 
GRH05A1 0 
GRH05A20 
GRH05A30 
GRH05AA0 
GRH05A50 
GRH05A60 
GRH05A70 
GRH05A80 
GRH05A90 
GRH05AA0 
GRH05AB0 
GRH05ACO 
GRH05AD0 
GRH05AE0 
GRH05AF0 
GRH05500 
GRH05510 
GRH05520 
GRH05530 
GRHO55A0 
GRH05550 
GRH05560 
GRH05570 
GRH05580 
GRH05590 
GRH055A0 
GRH055B0 
GRH055C0 
GRH055D0 
GRH055E0 
GRH055F0 
GRH05600 
GRH05610 
GRH05620 
GRH05630 
GRH056A0 
GRH05650 
GRH05660 
GRH05670 
GRH05680 


E-15 



c 

c 

c 


GO TO 510 


1700 N0SC#N0SC-1 _ 

IFSN0SC<1710, 1720, 1720 
1710 NOSC #0 
1720 CHISI<#CHIOLD 
1730 CONTINUE 


A SUCCESSFUL GIGANTIC STEP HAS OCCURRED. 

AN UNSUCCESSFUL GIANT STEP HAS OCCURRED. 
DELETE ITS OSCILLATION INFORMATION. 
NOSC#MAXOSNOSC-1,0< 


DECOMMENT FOR SSW OPTION 

CALL SSW TCH SNSSW,JUMP< 

I FS JUM P- 1 <2090 » 2090 ,1740 


# £ # # * * 


1740 IFSNF-NFMAX<1750 ,1750 ,207 0 
q 17^0 CONTINUE end qf the ma , n q0 LQOp> 

C ANOTHER CYCLE THROUGH THE VARIABLES HAS BEEN COMPLETED. 

C PRINT ANOTHER LINE OF TRACES. 

' 1760 jRfTE* ACE<1770 ’ 177 °wl 1AOO<CHIOLD,NF,%VEC2J<,J#1,NV< 

1770 I FSNZIP<18 10, 1780, 1810 
1780 IFSNTRACE <1810 , 1 810 , 1790 

1790 WRITES KW , 1410<%XSJ<, J#1 »NV< 

WRITES KW , 18 00< 

1800 F0RMATS1H < 

1810 NZIP#N7IP£1 
GO TO 520 

C A minimum has been found, print the remaining traces. 

C J5IS i F R?Ve T i ACE<1 ”‘ 0,1 ' 4 K 0 v!: S l2So<CM.tlLO,»F,,V E C lJ <. J ,l,.< 

WRITES KW,1410<SXSJ<,J#1,NV< 

C DECREASE THE SIZE OF THE STEPS FOR ALL VARIABLES. 

1840 CONTINUE DECOMMENT FOR SSW OPTION .... 

C CALL SSW TCH SNSSW , JUMP< 

C I F? JUMP- 1 <2090 , 2090 ,1 850 

^ 850 IFSNF-NFMAX<1860, 1860,2070^ A| _ L w%j< >le> delmin3;j< . 

I860 NGATEAl 

DO 1910 J# 1 , N V 

C if?mas ^ j<<1910 ' 1870 ’ 191 a°dx#abssdxsj« 

1870 ADX#DXSJ< 

I FSADX<1880, J890 ,1890 

1890 IFSADX-DELMINSJ<<19 10, 1910, 1900 

1900 NGATE#0 

1910 DXSJ<#DXSJ</RATIO 

c IFSNGATE<1920,1920,1960 check ^ JFLAJ%J<> 

1920 JFLM I N #5 

DO 1950 J#1,NV 
IFSMASKSJ <<1950,1 93 0,1 950 
1930 I FSJFLATSJ<— JFLMIN<1940 ,1950, 1950 
1940 JFLM IN #J FLAT SJ< 

1950 CONTINUE 

IFSJFLMIN-K2040 ,1990, 1960 
1960 IFSNTRACE<2130, 1970, 1970 

1980 F0RMATS///65H TERM* NATEDWHEN THE STEP SIZES BECAME AS SMALL AS 
*E DELM I N SJ < • < 

GO TO 2130 

1990 IFSNFLAT<2040 ,2040,2000 
2000 IFSNTRACE<2130, 2010, 2010 

2020 F0RMATS///72H TERM* nAtE DWIIEN THE FUNCTION VALUES AT ALL TRIAL 
♦NTS WERE IDENTICAL. < 

HRTTFS KW .2030<*DXSJ<, J«1,NV< _ 

2030 FORMATS/ / /2 9H CURRENT STEP SIZE VALUE S....//?1X9E13.5<< 

GO TO 2130 

c 2040 IFSNTRACE<510, 510, 2050 pRiNT ^ DXJJ< ^ search some more< 
C 2050 WRITES KW , 206 0<SDX% J< , J #1 , NV < 


GRH05690 
GRH056A0 
GRHO56B0 
GRH056C0 
GRH05600 
GRH056E0 
GRHO56F0 
GRH05700 
GRH0571 0 
GRH05720 
GRH05730 
GRH05740 
GRH05750 
GRH05760 
GRH05770 
GRH05780 
GRH05790 
GRH057A0 
♦ GRH057B0 
GRH057C0 
GRH057D0 
GRH057E0 
GRH057F0 
GRH05800 
GRH05810 
GRH05820 
GRH05830 
GRH05840 
GRH05850 
GRH05860 
GRH05870 
GRH05880 
GRH05890 
GRH058A0 
GRH058B0 
GRH058C0 
GRH058D0 
GRH058E0 
GRH058F0 
GRH05900 
GRH05910 
GRH05920 
GRH05930 
GRH05940 
GRH05950 
GRH05960 
GRH05970 
GRH05980 
GRH05990 
GRH059A0 
GRH059BO 
GRH059C0 
GRH059D0 
GRH059E0 
GRH059F0 
GRHO5A00 
GRH05A 1 0 
GRH05A20 
GRH05A30 
GRH05A40 
GRH05A50 
GRH05A60 
GRH05A70 
GRH05A80 
GRH05A90 
GRH05AA0 
GRH05AB0 
GRH05AC0 
GRH05ADO 
THGRH05AEO 
GRH05AFO 
GRH05B00 
GRH05B 1 0 
GRH05B20 
GRHO5B30 
POIGRH05B40 
GRH05B50 
GRH05B60 
GRH05B70 
GRH05B80 
GRH05B90 
GRH05BAO 
GRHO5BB0 
GRH05BC0 
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MORE THAN NFMAX # 17, 


2060 FORMAT*//60*1X1H*<//26H STEP SIZES REDUCED TO ....//% 1 X9 E 1 3 . 5<< 
GO TO 510 
C 

2070 WRITE* KW,2080<NFMAX 

2080 FORMAT *///46H ABNORMAL TERMINATION.... 

* 3 1 H CALLS TO THE CHIS0 SUBROUTINES 

GO TO 2110 

C DECOMMENT FOR SSW OPTION .... 

2090 WRITE* KW • 2 1 00< 

2100 FORMAT*///42H SUBROUTINE STEPIT TERMINATED BY OPERATORS 
C 

2110 WRITE* KW,2030<*DX*J<,J#1 ,NV< 


2120 KW I T #1 


SET SWITCH FOR TERMINATION. 

CALL FUNK WITH THE BEST SET OF X*J<-. 


2130 JVARY #0 

CALL FUNK 

I F*CHISQ-CHISAV<2150,2150,2 140 
2 140 WRITE *KW ,410< CH ISA V ,CH I SO ,NF 
2150 IF*NTRACE<2200, 2160, 2160 
2160 WRITE* KW ,217 0<NF 

2170 FORMAT*// 1X,I5,23H FUNCTION COMPUTATIONS < 

WRITE* KW ,2180<*X*J<, J# 1 ,NV< 

2180 FORMAT*/ // 10X2 4HF I NAL VALUES OF X * I <....//* 1 X 5E 2 1 . 1 3« 
WRITE* KW .2190<CHIS0 

2190 FORMAT*/ / 24H FINAL VALUE OF CH I SO # E21.13//< 

2200 IF*KWIT<2250, 2210, 2250 

MATD#IABS*MATRIX-100< 

2210 MATO AMATRIX— 100 

I F*MATD<2220, 2230, 2230 
2220 MATO#— MATO 

2230 I F*MATD-50<2240, 2240, 2250 


SKIP ERROR CALCULATION 
2240 IF*NACTIV-NV<2250, 2260, 2250 
2250 RETURN 


IF ANY MASK*JC.NE.O 


2260 CONTINUE 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

TO DELETE THE ERROR COMPUTATION SECTION, REMOVE ALL STATEMENTS 
BETWEEN HERE AND THE END STATEMENT, AND INSERT A RETURN STATEMENT. 

COMPUTE THE STANDARD ERRORS AND THE CORRELATIONS. 

FIRST COMPUTE THE SECONO PARTIAL DERIVATIVES OF CH I SO WITH RESPFCT 
TO THE X*J<. 


FAC#RATI0***MATRIX-100< 
00 2300 I #1 ,NV 


COMPUTE THE DIAGONAL PARTIALS FIRST. 
DX*I <#ABS*FAC*DX*I« 


DX*I<#FAC*DX*I< 

I F*DX*I <<2270, 2280, 2280 
2270 DX* I <#-DX *1 < 

2280 X SA VE *1 <#X *1 < 

JVARY #0 

DO 2290 J# 1 , 2 

X*I<#XSAVE*I<6DX*I< 

CALL FUNK 
NFANFC 1 
J VAR Y # I 

SEC0ND*1 , J<#CHI SO 
2290 DX*I <#— DX*I < 

X* I <#X S A VE *1 < 

2300 ERR* I ,I<#*SEC0ND*1 , 1 <-2 . 0*CH I OLD&SECOND* 1 ,2«/DX*I<**2 


IF*NV-2<2350 ,2310,2310 
2310 DO 2340 I#2,NV 
I M # I — 1 

DO 2340 J# 1 , 1 M 

DO 2330 K#1 ,2 

X*I <#X SA VE*I <£DX*I < 

JVARY #0 

DO 2320 L# 1 , 2 

X* J <#X S A VE*J <£DX*J< 

CALL FUNK 


COMPUTE THE OFF-DIAGONAL PARTIALS. USE A 
REDUNDANT FOUR-POINT RULE FOR GREATER 
RELIABILITY. 


GRH05BD0 
GRHO5BE0 
GRH05BF0 
GRH05COO 
GRH05C10 
GRH05C20 
GRH05C30 
GRH05C40 
GRH05C 50 
GRH05C60 
GRH05C70 
GRH05C80 
GRH05C 90 
GRH05CA0 
GRH05CB0 
GRH05CC0 
GRH05CD0 
GRH05CEO 
GRH05CF0 
GRH05000 
GRH05D10 
GRH05D20 
GRH05D30 
GRH05D40 
GRH05D50 
GRH05D60 
GRH05D70 
GRH05D80 
GRH05D90 
GRH05DA0 
GRH05DBO 
GRH05DC0 
GRH05DD0 
GRH05DE0 
GRH05DFO 
•GRH05E00 
GRH05E 10 
GRH05E20 
GRH05E30 
GRH05E40 
GRH05E50 
GRH05E60 
GRH05E70 
GRH05E80 
GRH05E90 
GRH05EA0 
GRH05EB0 
GRH05 ECO 
GRH05EDO 
GRH05EE0 
GRH05EF0 
GRH05F00 
GRH05F1 0 
GRH05F20 
GRH05F30 
GRH05F40 
GRH05F50 
GRH05F60 
GRH05F70 
GRH05F80 
GRH05F90 
GRH05FA0 
GRH05FB0 
GRH05FC0 
GRH05FD0 
GRH05FE0 
GRH05FF0 
GRH06000 
GRH06010 
GRH06020 
GRH06030 
GRH06040 
GRH060 50 
GRH06060 
GRH06070 
GRH06080 
GRH06090 
GRH060A0 
GRH060BO 
GRH060C0 
GRH060D0 
GRH060E0 
GRH060FO 
GRH06100 
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NF #NF6 1 GRH06110 

J VAR Y # J GRH06120 

S ECONO *K,L<#CHISQ 959n5wn 

X*J<#XSAVE*J< GR HO 61 AO 

2320 DX*J<#-DX*J< GRH06150 

X*I <#XSAVF*I < GRH06160 

DECOMMENT FOR SSW OPTION .... GRH06170 

CALL SSWTCH %NSSW,JUMP< 55925 , 1 52 

1 F* JUMP- 1 <2090, 2090, 2 330 GRH06190 

nx*T<#-DXXI< GRH061B0 

- ERR* I tJ<#0.25**SEC0ND*l, 1C-SEC0ND* 1,2<-SEC0ND*2 , 1<G SECONDS2 , 2« GRH061C0 

* /*DX*I <*DX*J« GRHO61D0 

?740 FRR % I . T < AF RR X I • J < GRHO61E0 

2340 tWRXJ f Q(= 0ERIVATIVE COMPUTATION. GRH061F0 

GRH06200 

2350 IF*NTRACE<2410, 2360, 2360 ' 95995^9 

WRITER KW*2370< GRH06220 

2370 FORMAT *4 1H1SIZES OF INCREMENTS TO BE USED BELOW.. ..< GRH06230 

WRITE* KW,2380<*DX*J<, J#1 ,NV< 25995?£2 

2380 F0RMAT*/*1X9E 13 .5« GRH06250 

WRITE* KW , 2 39 0< GRH06260 

2390 F0RMAT*/////45H MATRIX OF THE SECOND PARTIAL DERIVATIVES.... /< GRH06270 

DO 2400 I # 1 , N V 25995?£2 

2400 WRITE* KW,2380<*ERR*I ,J<,J#1,I< GRH06290 

2410 DO 2420 I#1,NV GRH062A0 

DO 2420 J#1,I 95995?r9 

I F*ERR*I ,J«2420, 2430, 2420 25.995I92 

2420 CONTINUE GRH062D0 

GO TO 2450 GRH062E0 

2430 WRITE* KW,2440< GRH062F0 

2440 F0RMAT*////46H THE ABOVE MATRIX CONTAINS ONE OR MORE ZEROES. / GRH06300 

* 75H A LARGER VALUE OF -MATRIX- SHOULD BE TRIED, TO SEE IF THEYGRHO6310 

* ARE LEGITIMATE. < 25995212 


* ARE LEGITIMATE. < GRH06320 

GRH06330 

j,**#####*##*:*:****#*****’! 1 *******#*** GRH06340 

GRH063 50 

INVERT THE MATRIX OF SECOND DERIVATIVES USING THE ALGORITHM SYMINV2 GRH06360 
*ALG0R I THM 150, H. RUTISHAUSER, COMMUNICATIONS OF THE A.C.M. 6 * 1 963<GRH063 70 
P. 67< . RE-USE SOME STORAGE. GRH06390 

2450 SGNDE T # 1 • 25995?£2 

DETLOGSO. 95995*59 

DO 2460 J# 1 »NV 25995222 

2460 SAL VO*J <#1 .0 95995??9 

nn 2640 T*1, M V GRH063E0 

’ " PFRFORM PIVOT SEARCH. GRH063F0 

BIGAJJ#0.0 25995/ 4 92 

DO 2510 J#1,NV GRH06410 

IF*SALV0*J«2470, 2510, 2470 „ 25995512 

ABER#ABS*ERR*J»J<< GRH06430 

2470 ABER #ERR * J » J< 25995552 

IF*ABER<2480 ,2490 ,2490 GRH06450 

2480 ABERS-ABER , 2599 5552 

2490 IF*ABER-BIGAJJ<2510, 2510, 2500 GRH06470 

2500 BIGAJJKABER GRH06480 

K#J GRH06490 

2510 CONTINUE 259955£2 

IF*BIGAJJ<2540,2520,2540 25995599 

2520 WRITE* KW,2530< GRH064C0 


2500 BIGAJJKABER 259955?2 

K#J GRH06490 

2510 CONTINUE 259955£2 

IF*BIGAJJ<2540,2520,2540 25995599 

2520 WRITE* KW,2530< GRH064C0 

2530 FORMAT *////67H ERROR MATRIX IS SINGULAR. -MATRIX- SHOULD PROBABL YGRH064D0 

* BE INCREASED. /////< 25995512 

GO TO 2120 GRH064F0 

2540 SALVO*K<#0 • 0 25995192 

I F*E RR *K ,K<<2 550,2 55 0,256 0 GRH06510 

2550 SGNDE T#- SGNDE T GRH06520 

C2 560 DETLOG#DE TL0G6AL0G*B IG AJJ </2 .303 25995119 

2560 DETLOG#DETLOG6DLOG*BIGAJJ</2.303 25995152 

r * GRH06550 

c PREPARE FOR THE NEXT ELIMINATION STEP. GRH06560 

TR I AL*K<#1 ,0/ERR*K,K< 95999519 

ERR*K,K<#0 . 0 25995122 

XS A VE *K <#1 .0 25995112 

m#k-1 GRH065A0 

IF*M<2600, 2600, 2570 259n5trn 

2570 DO 2590 J#1,M GRH065C0 

xsa ve *j <#e rr *k , j< 95995199 

TRIAL *J<#ERR*K,J<*TRIAL*K< 25995512 

1F*SALV0*J <<2520, 2590, 2580 95995559 

2580 TRI AL*J<#-TRIAL*J< 95995592 

2590 ERR*K , J <#0 .0 259955*9 

2600 M#KG1 GRH06620 

I F*M-N V<26 10 ,2610,2650 95995539 

2610 DO 2640 J#H,NV GRH06640 
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XSAVESJ<#ERRSJ,K< „ „ 

I ESS AL VOS J<<2 520 ,2620,2630 

2630 TRI ALSJ<#-ERRSjfK<*TRI ALSK< 

^2640 ERRSJ,K<#0.0 PERFORM THE NEXT ELIMINATION STEP. 

2650 DO 2660 J#1,NV 
DO 2660 K#J,NV 

2660 ERRSK, J<#ERRSK , J <£X S A VE SJ < *TR I ALSK< 

PRINT THE ERRORS AND CORRELATIONS, AND TERMINATE. 


IS NEGATIVE DEFINITE. 


IFSSGNDET<2670 ,2690 ,2690 
2670 WRITES KW,2680< 

2680 FORMATS////75H ERROR MATRIX 
* PROBABLY BE DECREASED. < 

2690 IFSNTRACE<2720, 2700, 2700 

2710 F0RMATS////3 9H < L0G10 F S^ETERM I N AN T OF ABOVE MA TR I X< # E12.5.10X 


GRH06650 
GRH06660 
GRH06670 
GRH06680 
GRH06690 
GRH066A0 
GRH066B0 
GRH066CO 
GRH066D0 
GRH066E0 
GRH066F0 
GRH06700 
GRH06710 
GRH06720 
GRH06730 
GRH06740 
-MATRIX- SHOULDGRHO6750 
GRH06 760 


- - 22HSIGN OF DETERMINANT E * E M< ^ TWICE TH g INVERSE 0 F 

; THE MATRIX OF SECOND DERIVATIVES. 

2120 DO 2790 I #1 ,NV 
DO 2730 J#1,I 
ERRSI , J <#E RR S I , J<*2 .0 
2730 ERRSJ,I<#ERRSI ,J< 


GRH06770 
GRH06780 
GRH06790 
GRH067A0 
GRHO67B0 
GRH067C0 
GRHO67D0 
GRH067E0 
GRH067F0 
GRH06800 

XSAVESI<#SIGNSSQRTSABSSERRSI ,I<«, ERRSI, I<<GRH06810 
BtxatKKsi ,iv GRH06820 

I FSABER<2 740 ,2790,2750 
2740 ABER#-ABER 
C2750 ABER#SQRT SABER< 

2750 ABER#DSQRTSABER< 

I FSERR SI , I <<2760,2770,2790 
2760 ABER#-ABER „ 

2770 WRITES KW , 278 0<E RR SI , I < 

2780 FORMATS ///50H NEGATIVE OR ZERO MEAN SQUARE ERROR ENCOUNTERED... 

* 3XE15.8/39H -MATRIX- SHOULD PROBABLY BE DECREASED. ///< 

2790 X SAVE SI <#ABER 

I FSNT RACE <21 30 ,2800 ,2800 

2800 WRITES KW,2810< 

2810 F0RMATS/////20H STANDARD ERRORS.... < 

WRITES KW , 2 38 0<SX S A VES J < , J # 1 ,NV< 

I FSNV— 1<2120, 2120,2820 
WRITES KW , 2830< 

F0RMATS/////45H LOWER TRIANGLE OF THE CORRELATION MATRIX. 

DO 2880 I #2 , N V 
I M # I — 1 

DO 2870 J#1,IM 
DENOM #XSA VE SI <*XSAVESJ< 

IFSDENOM<2850, 2840, 2860 
2840 TRIAL SJ<#0. 

GO TO 2870 
2850 DENOM #-DENOM 
2860 TRIALSJ<#ERRSI ,J</DENOM 

2870 CONTINUE , 

2880 WRITES KW , 2 38 0<STR I AL S J<, J # 1 , I M< 

GO TO 2120 
END 

BLOCK DATA 


2820 
2 830 


/< 


BLOCK DATA SUBPROGRAM FOR STEPIT, SIMPLEX, 
J. P. CHANDLER, F.S.U. PHYSICS DEPT. 


AND STP. 


/* 

// 


COMMON /FRODO/ NFM AX ,NFL AT , J VAR Y ,NEXTR A 
COMMON/ HISTO / I CUR VE ,CR VM AX , A BE S 1 02< , BI NS 1 02< 
DATA NFMAX/1000000/ , NFL AT/ 1/ , NEXTRA/O/ 

DATA ICURVE/O/, CRVMAX/O./ 

END 


GRH06830 
GRH06840 
GRH06850 
GRH06860 
GRH06870 
GRH06880 
GRH06890 
GRH068AO 
GRHO68B0 
GRH06BC0 
GRHO68D0 
GRH068EO 
GRH068F0 
GRH06900 
GRH06910 
GRH06920 
GRH06930 
GRH06940 
GRH06950 
GRH06960 
GRH06970 
GRH06980 
GRH06990 
GRH069A0 
GRH069BO 
GRH069C0 
GRH069D0 
GRH069E0 
GRH069F0 
GRH06A00 
GRH06A10 
GRH06A20 
GRH06A30 
GRH06A40 
GRH06A 50 
GRH06A60 
GRH06A70 
GRH06A80 
GRH06A90 
GRH06AA0 
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E. 5 SUBROUTINE USED FOR GENERATING A GREAT CIRCLE ARC 


MAIN 

THIS PROGRAM GENERATES A TEST ARC OF VARIABLE LENGTH IN 
CARTESIAN COORDINATES, PROJECTS AND TRANSFORMS TO DATA BASE 
COORDINATES AND COMPARES DIREC1 AND INTERPOLATED VALUES FOR 
INTERMEDIATE POINTS ON THE ARC. 

PROGRAMMER - E. MICHAEL O'NEILL 

COMPUTER SCIENCES CORPORATION 

8728 COLESVILLE RD. SILVER SPRING, MD. 


DIMENSION X Y ZO ( 3 ) , XYZ(3), BXYZI3) 

DATA XYZO ,P 1/1. 0,0.0 , 0.0,3. 14159/ 

TARC=0.02 
DARC=P 1/16.0 
ANG=0 . 0 
ARC=0 .0 
BANG=P 1/2.0 
DO 10 1=1,5 

CALL CIRARC(XYZO,ANG»ARC»BXYZ) 

ARK = 0 . 0 
DO 5 J= 1 » 5 
WRITE(6,4) I,J 

4 FORMAT ( ' I , J ' ,2 13/ ) 

CALL CIRARCIBXYZ, BANG, ARK, XYZ) 

CALL TSTARC(XYZ,TARC) 

5 ARK= ARK+DARC 
10 ARC=ARC+DARC 

STOP 

SUBROUTINE TSTARC (XYZO,ARCO) 

DIMENSION XYZO (3 ) , XYZ1 ( 3 ) ,XYZ ( 3) 

DATA PI/3.14159/ 

DANG=P 1/16.0 
0ARC=ARC0/4.0 
ANG=P I 

DO 20 I ANG= 1 , 3 

CALL CIRARC ( XYZO.ANG, ARC0,XYZ 1) 

CALL QTRAN(XYZO,PO,EO, IFO) 

CALL QTRAN(XYZ1,P1 , E 1 , I F 1 > 

DP=( PI— PO 1/4.0 
DE=( E1-E0) /4.0 
P=PO 
E = EO 
A RC = 0 . 0 

DO 10 I ARC= 1,5 

CALL CIRARC (XYZO, ANG, ARC, XYZ) 

CALL OT RAN (XYZ, PC, EC, IFC) 

WRITE (6,5 i PC » P »FC,F , IFC, XYZ 
5 FORMAT! 1X,4F12.3, I 12,3F12.3) 

P= P + DP 
E=E+DF 

10 ARC= ARC+DARC 
WRITE! 6,5) 

20 ANG= AN G+DANG 

WRITF(6,5) 

RETURN 

END 

SUBROUTINE C I R A RC ( X Y ZO , ANG , ARC , X YZ ) 

DIMENSION XYZO ( 3 ) , XYZ ( 3 ) 

RAD=SORT ( XYZO ( 1 )**2+XYZ0 ( 2 ) **2+XYZ0 ( 3) **2) 

RHO = S OR T ( XYZO ( 1 ) **2+XYZ0 ( 2 ) **2 ) 

C1=C0S( ANG) 

S 1 = S I N ( ANG) 

C2=C0S( ARC ) 

xYZUUixy/n ( i )4 xyzd ( 3 )*c i-r ad*xyzo ( 2 ) *si )*s2/rho + 

XYZ ( 2 )= (XYZO ( 2)*XYZ0 ( 3 ) *C1 + RAD*XY Z0( 1 ) *S1 )*S2/RH0 + 

X Y Z ( 3 > = -RH0*C1*S2 + XYZO ( 3 ) *C2 

RETURN 

END 


XY Z 0 ( I ) * C 2 
X Y ZO ( 2)*C2 
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E, 6 SUBROUTINES USED FOR MAPPING A POINT IN CARTESIAN SPACE 
INTO DATA BASE COORDINATES 


SUB ROUT IMF OTRAN (COORD ,X , Y , I FACE ) 

THIS SUBROIIT IMF TAKES AS INPUT THE CARTESIAN COORDINATE ARRA-Y 
(COORD) AMD RETURNS DATA BASF COORDINATES (X AMD Y) AND THF 
FACE NUMBER ( I FACE ) OF THE X,Y SYSTEM 

PROGRAMMER - E. MICHAEL 0 ' ME I L L 

COMPUTER SCIENCES CORPORATION 
8728 COLES V I LLE RD. SILVER SPRING, MO. 

COMMON /TRACOM/ RS , R 0 , R S O , R 4 , G AMM A , OM G AM , l)E L T A , OM EG A , 

1 COO , C 1 0 , CO 1 , C 1 1 , C 20 , CO 2 , DO , D 1 , 02 

DIMENSION CnORDm, C<3), ROT ( 3 ,/> ) , IROT(3,6) 

DATA ROT/1. 0,1. 0*1.0, - 1 .0 , 1 . 0 , - 1 . 0, - 1 . 0 , 1 . 0 , 1 . 0 , 

1 1 . 0 , 1 . 0 , - 1 . 0 , -1.0, -1.0, 1.0, -1.0, 1.0, -1.0/ 

DATA TROT /?., 3, 1,2, 3, 1,1, 3, 2,1 ,3, 2, 1,2, 3, 1,2,3/ 

DATA OFFSET/2048. 0/ 

C DETERMINE FACE NUMBER BY FINDING LARGEST (ABS) COMPONENT 

CMAX=ARS(COORD ( 1 ) ) 

J=1 

on 10 1=2,3 

IFICMAX.LT. (ABS(CO(IRD( I ))) ) J=I 
CMAX=ARS(COORD ( J ) ) 

10 CONTINUE 
IFACE=J*2 

I F ( COORD ( J ) . GE . 0 . 0 ) I FACE = I FACE-1 
C INTERCHANGE CARTESIAN COMPONENTS TO GET TO FACE SYSTEM 

DO 20 1=1,3 

20 C ( I ) = C OOR 0 ( I ROT ( I , I FACE ) ) *ROT ( I , I FACE ) 

C PROJECT TO PSI, F 1 A SYSTEM 

PRJ=R0/C(3) 

P = C ( 1 ) *P R J 
E=C(2)*PRJ 

C TRANSFORM TO X,Y 

X= FM A P ( P ,E ) /RO*OFFSFT + OFFSET 
Y = FMAP(E,P) / R0 s;! 0 F F S E T + 0 F F S £ T 
RETURN 
END 
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FUNCTION FMAP ( A ,B ) 

THIS FUNCTION PERFORMS THE INVERSE TRANSFORMATION (F STAR) 


C 

c 


PROGRAMMER - E. MICHAEL O'NEILL 
K ’ COMPUTFR SCIENCES CORPORATION 

872 B COLESVILLE RD. SILVER SPRING, MD. 

COMMON /TRACOM/ RS , RO , RSO , R4 , GAMMA, OMGAM, DELTA , OMEGA, 

1 COO, Cl 0, CO 1,C11,C20,C02, 00,01,02 

ASO= A*A 
A3=A*AS0 
A4= ASO*A SO 
BSO=B*B 
B4=BSO*BSO 

FMAP= A#GAMMA+A3# ( OM GAM/RSO ) +A*BSQ*R MA* ( D EL T A+ ( R SQ-B SO ) * 

1 ( C00+C10*ASQ+C01*BSO+C 11*ASQ*BSQ+C20*A4+C02*B4) ) + 

2 A3*RMA* (OMEGA+RM A* (DO+D 1 * A S0+ 02* A4 ) ) 

RETURN 

ENO 

B DATA INITIALIZATION FOR FMAP 

COMMON /TRACOM/ RS , R 0 , R SO , R4 , G AMM A , CJMGAM , DEL T A , OM EGA , 

1 C00,C10,C01 ,C11 ,C 20, CO 2, DO, 01 ,D2 

S*"‘ cS8,f{0,C0i;cll,C 2 0.C0 2 M.»00 ? || ; _ 2 |.47 ?U j-«.9|«^ 
DATA 00 , 01 , 02/-I2. 80 531, 50. A 65 1 1,-256.5819/ 
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E.7 SUBROUTINES USED FOR SIMULATING THE FAST-FILLING METHOD 


SUBROUTINE SCAN1(INBASE,NSL,NBL) 


THIS SUBROUTINE SIMULATES THE STORAGE OF DATA IN THE DATA BASE 
DM A SINGLE INTERPOLATION BLOCK BASIS. 

INTERPOLATION ALONG AND BETWEEN SCAN LINES TAKES PLACE USING 
VFCTOR COMPONENTS INPUTTED ( I P Z , I E Z , I D P . I DE , JDP , J DE ) . Akir> 

THE SERIAL ADDRESS IS CALCULATED USING A TABLE LOOK UP AND 

THE RECORD NUMBER IS DETERMINED BY A SHIFT OF THE SERIAL ADDRESS, 


PROGRAMMER 


E. M ICHAEL O'NEILL 

COMPUTER SCIENCES CORPORATION . 

8728 COLESVILLE RD . SILVER SPRING, MD. 


C 


// 

// 


LOGICAL FAILW 
LOGICAL*]. INDATA, IODATA 

INTEGER SHETL , SHFTR , 

COMMOM/DATCOM/ INDATA ( 81 920) , I ODAT A ( 24 5 76) 
COMMON / I NDCOM / INDX<4096), INDY (4096), ISHFT, 
COMMON /SC N COM/ IPZ, IEZ, IOP, IDE, JDP, JDE , 
DATA LREC/-1/ 

KLOC= INBASE 
DO 20 I B L = 1 , N B L 
JLOC=KLOC 
IP= IPZ 
I E = I E Z 

DO 10 I S L = 1 , M S L 

T = T P / ISHFT 

J = JP/ ISHFT 

KFY= I M D X ( I ) + I ND Y ( J ) 

IRFC = SHFTR (KEY, 6) 

I F ( IRFC.FO.LREC ) GO TO 5 
LR FC= I RE C 

CALL I OMMGR ( IREC, IOBA SE , F A I LW ) 

5 I F ( FAILW) GO TO 9 

I LOC = KFY-SHFTL ( IRFC,6) 

IDUM= IOBASE + IL0C 

I OOA T A ( IOBASF+ I LOC ) = I NDA TA ( JLOC ) 

9 IP=IP+IDP 
10 I E = I F+ I DF 

KLOC=KLOC+ JSIZE 
IP 7= IPZ + JOP 
20 I F Z = IFZ + JDF 
RETURN 

FNn 

SUB ROUT INF IOMNGR ( I R EC , 1 0 B A SE , F A I LW ) 

LOGICAL FAILW 

IOBASF=0 

FA ILW= .FALSE . 

RF TURN 
END 

FXFC. LOADER , REGION = 2 OOK ,PARM = ' SI ZE = 200000' 


NMI L 
JSIZE 
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E. 8 SUBROUTINES USED FOR TESTING I/O HANDLING METHOD 


SUBROUTINE IDMNGR ( IKEY, IOBASE»FAILW) 

THIS SUBROUTINE CONTROLS DATA BASE INPUT/OUTPUT 


THIS SUBROUTINE TAKES AS INPUT THE SERIAL ADDRESS OF A 
RECORD REQUIRED FOR PROCESSING, COMPARES THF RECORD NUMBFR 
WITH A LIST OF RECORDS AVAILABLE IN CORE BUFFER AND FFTCHFS 
NFW RECORDS WHEN NEEDED. 

THE CORE BUFFER IS CONTROLED BY A PUSH DOWN STACK, 

WHEN A NEW RECORD IS REQUIRED 1 HE OLDEST RECORD IN THE BUFFFR 
IS RETURNED TO THE DATA BASE AND REPLACED BY THE REQUESTED " 
RECORD. THE BASE ADDRESS OF THE REQUESTED RECORD IS RETURNFD 
TO THE CALLING PROGRAM IN IOBASE. 

IOMNGR ALSO RETURNS THE WINDOW FLAG TO INDICATE PRESENCE 
OR ABSENCE OF THE REQUESTED RECORD IN THE DATA BASE. 


PROGRAMMER - E. MICHAEL O'NEILL 

COMPUTER SCIENCES CORPORATION 

8728 COLESVILLE RD . SILVER SPRING, MD. 

LOGICAL FAILW, FAIL(3,6) 

DIMENSION INDEX (3,6) , LNDEX ( 3 ) 

EQUIVALENCE (FAIL! 1,1) , INDEX! 1,1) ) 

DATA NSTAK/O/, MSTAK/6/ 

DATA INDEX /— 1,6* 0,-1, 5,0, — 1,4,0,— 1,3,0, - 1,2,0,- 1,1,0/ 
ON FIRST ENTRY IMMEDIATELY GFT REQUESTED RECORD 


IF! NSTAK.EO.O) GO ID 25 

TEST REQUEST RECORD NUMBFR AGAINST RECORDS IN CORE 

DO 20 I=1,NSTAK 
JKFY= INDEX ( 2,1) 

BRANCH IF A MATCH IS FOUND 

IF! IKEY.EO. INDEX! 1, I ) ) GO TO 30 
20 CONTINUE 

BRANCH IF STACK IS NOT YET FULL 
IF ( NSTAK.LT. MSTAK) GO TO 25 


ENTER REQUESTED RECORD INTO STACK 


CALL PUSHD ( INDEX, LNDEX, MSTAK, IKEY) 
JKFY= LNDEX ( 2 ) 


RETURN OLDEST RECORD TO DATA RASE 
CALL PUT RFC (LNDEX ( 1 ) , LNDEX (2 ) ,F AILW) 


RFAD MEW RECORD INTO SPACE VACATED 

CALL GET RFC ( INDEX ( 1 , 1 ) , INDEX ( 2, 1 ) , F A I LW ) 

FAIL(3,JKEY)=FAILW 

GO TO 30 

FNTFR AND FFTCH ONLY UNTIL STACK IS FULL 


25 N S T A K = N S T A K + 1 
•JKEY = NSTAK 

CALL PUSHD! INDE X , LNDE X , MS T A K , I KE Y ) 

CALL GETREC ( IKEY , JKEY,F A JLW ) 

FA IL ( 3 , JKEY ) = F A I L W 

SET BASE ADDRESS AMD WINDOW FAILURE FLAG BEFORE RETURNING 

30 IORASE= ( JKFY-1 ) *1024 
FA ILW = FA IL (3 , JKEY) 

RETURN 

END 
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SUBROUTINE PUSHD! I S T A CK , L , M AX , I IM ) 


PUSH DOWN STACK ROUTINE 


CONTROLS TRIPLE STACK CONTAINING CURRENTLY AVAILABLE RECORD 
NUMBERS AND THE ASSOCIATED BUFFER POINTER AMD WINDOW FLAG 
ISTACK IS AN ARRAY CONTAINING THE STACK U 

L IS THE VECTOR ELEMENT OF ISTACK BEING REPLACED BY IN 
MAX IS THE DEPTH OF THF STACK 
IN IS THE NEW SCALAR ADDED TO THE STACK 
IN IS AUTOMATICALLY ASSOCIATED WITH THE POINTER FREED BY L 


PROGRAMMER - E. MICHAEL O'NFTLL 

COMPUTER SCIENCES CORPORATION 
8 72 R COLESVILLE RD. SILVER SPRING, MD. 

DIMENSION TSTACK(3,1), L ( 3 ) 

L ( 1 ) = I ST ACK ( 1 , M A X ) 

L ( 2 ) = I S T A CK (2, MAX) 

L ( 3 ) = I ST ACK ( 3 , MAX ) 

K=MAX- I 
J = M A X 

DO 10 1 = 1, K 

I STACK ( 1 , J ) = ISTACK ( 1 , J-l ) 

ISTACK! 2»J) = ISTACK(2»J-1) 

ISTACK! 3, J ) = I S TACK (3, J-l! 

10 J=J-1 

ISTACK! 1 , 1 ) = IN 
ISTACK! 2, 1 ) =L (2 ) 

RETURN 

END 

SUBROUTINE GET RFC (NREC, IPOS,EAILW) 


THIS SUBROUTINE HANDLES 1 HE ACTUAL READING AND WRITING OE 
DATA BASE RFCORDS. NREC, THF RECORD NUMBER AND IPOS, THE 
BUFFER RASE ADDRESS A R F INPUT. SUBROUTINE WINDOW IS CALLFD TO 
ASCERTAIN WHETHER OR NOT THE RECORD IS IN THF DATA BASE AMD ITS 
ABSOLUTE ADDRESS, MREC. FAILW, THE WINDOW FLAG, IS RETURNFD 
TO THF CALLING PROGRAM. 
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LOGICAL FAILW 
COMMON I DATA (6 14 A ) 

DATA N DBASE / 11/ 

CALL WINDOW TO OBTAIN ABSOLUTE RECORD NUMBER AND WINDOW FLAG 

CALL WINDOW (NREC, MRFC, FAILW) 

RFTURN IF RECORD DOES NOT FXIST 
IF! FAILW) RETURN 

CALCULATE BASE ADDRESS AND RFAD RECORD 
I = ( I PO S- 1 ) * 1024 + 1 
L = IPOS +1023 

RFAD! MORA SF 'MREC ) ( T D ATA ( J ) , J = J , L ) 

RFTURN 

ENTRY PUT RFC (MREC , IPDS ,FA I LW ) 

CALL WINDOW TO OBTAIN ABSOLUTE RFCOR.D NUMBER AND WINDOW FLAG 

CALL W TNDOW (MREC , MRFC , FA TLW ) 

C RFTURN IF RECORD DOES MOT FXIST 

IF! FA TLW) RFTURN 

C CALCULATE BASE ADDRESS AND WRITE DIJT RECORD 

I = ( I PO S- 1 ) * 10 24 + 1 
L= TP0S+1023 

WRITE! NDBASE 'MREC ) ( ID ATA ( J ) , J = I , L ) 

RETURN 
F MD 


k. 
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SUB RDIIT I ME WINDOW (NREC,MREC,FAILW) 

THIS SUBRflUT I NE WINDOWS THE DATA BASE AND ASSIGNS INDIRECT 
RECORD NUMBERS TO THE ACTUAL PHYSICAL RECORDS 

A CALL TO SETWIN INITIALIZES I HE RECORD NUMBER (IREC) AMD 
FAILURE FLAG (WIND) TABLES. 

A CALL TO WINDOW WITH NREC AS INPUT RESULTS IN THE RETURN 
AM ABSOLUTF RECORD NUMBER (MRFC) AND A WINDOW FAILURE FLAG 
( FA ILW) . 
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LOGICAL FAILW 
LOGICAL *1 WIND! 1536) 

DIMENS ION IREC ( 1536) 

DATA WIND/1536*. TRUE./ 

DATA I RE C / 1 53 b*0/ 

FA ILW=. FALSE . 

C S FT FLAG TO INDICATE RFC ODD IS PRESENT 

r. CAICHIATF WHICH RFCORD ri.ock rfodested record belongs to 

J R E C = ( NREC- 1 ) /16+1 
KRFC=NPEC-JRFC*16 

C LOOK UP ADDRESS IN ABSOLUTF ADDRESS TABLE 

MRFC=IREC( JRFC)*16+KREC 

C SET FAILURE FLAG IF RECORD DOFS NOT EXIST 

I F ( W IND ( JRFC ) ) FA I LW= .TRUE. 

RETURN 

ENTRY SETWIN 

C SET UP DATA BASF AREAS TO BE STORED, WIND=. FALSE. MEANS 

C THAT THE RECORD BLOCK EXISTS. 

WIND! 1 )=. FALSE. 

WIND! 2 ) = . FALSE. 

WIND! 5 )=. FALSE. 

WIND! 6 ) = . FALSE. 

C DERIVE ABSOLUTE ADDRESS TABLF FROM TABLE OF MON EXISTANT RFCORDS . 

J = 0 

DO 10 1=1,1536 
IF ( .NDT .WIND (I) ) J = J+ 1 
IF ( .MOT .WIND ( I ) ) IREC ( I )=J 
10 CONTINUE 
RETURN 
END 
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E. 9 SUBROUTINES USED FOR RETRIEVING DATA FROM THE DATA BASE 
AND FOR GRAPHIC DISPLAY 


SUBROUTINE R A S TER ( I SE R , I BLOK ) 

THIS SUBROUTINE UNPACKS FROM SERIAL TO RASTER FOR THE SIXTH 
LEVEL OF DIVISION ( 4096 - 64 X 64 ) 


THE METHOD EMPLOYED IS TO LOOP OVER THE DESIRED NUMBER OF 
LFVELS ADDING AT EACH LEVEL AN X AND A Y DISPLACEMENT 
APPROPRIATE TO THAT LEVEL AND OLJADRANT. 
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LOGICAL *1 I S E R ( 40 96 ) , IBL0K<64,64) 

DIMENSION 11(4), 12(4), 13(4), 14(4), 15(4), 16(4 
DIMENSION J1 (4) , J 2 ( 4 ) , J3 (4 ) ,J4(4),J5(4),J6(4) 

DATA 11,12,13,14,15,16 „ , , , „ , , „ , . . 

1/1,— 1,1, 0,1, — 3, 1,0,1, — 7, 1,0,1, 15,1,0,1, 31, 1,0, 1 , 63,1,0/ 

l/oIl,ol6^-{^lJ-l’ 0r3O,-3, 0,-7, 1,-7, 0,-15, 1,-15, 0,-31, 1,-31, 0/ 
1 = 1 
J=1 
N= 1 

LOOP OVER ALL LEVELS 


DO 60 L6= 1 ,4 
DO 50 L 5= 1 , 4 
DO 40 L4 = 1 , 4 
DO 30 L 3= 1 , 4 
DO 20 L2 = l ,4 

STORE ITEM FROM FIRST QUADRANT 

I B L OK ( I , J ) = I SER ( N ) 

N=N+ 1 


STORE ITEMS FROM QUADRANTS 2-4 

DO 10 L 1= 1 , 3 
I=I+I1( LI) 

J=J+J1( LI ) 

I BLOK ( I , J) = ISER (N ) 

10 N = N + 1 

I = I + I 2 ( L 2 ) 

20 J= J+J2( L2 ) 

I = I + 1 3 ( L3 ) 

30 J = J+ J3( L 3 ) 

I = I + I 4 ( L4 ) 

40 J= J+J4( L4 ) 

I = I + I 5 ( L 5 ) 

5Q J= J+J5( L5 ) 

I = I + I 6 ( L6) 

60 J= J + J6( L 6 ) 

RETURN 

END 
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SUBROUTINE GRA Y3 ( NGRA , I 1 28 ) 


THIS SUBROUTINE CREATES A GRAY SCALE PLOT OF A 128 X 128 
INPUT ARRAY ON THE PRINTER. GRAY 3 OVERPRINTS UP TO TWO TIMES 
FOR TONES RANGING FROM WHITE (1) TO BLACK (16). 

PROGRAMMER - E. MICHAEL O'NEILL 

COMPUTER SCIENCES CORPORATION 

8728 COLESVILLE RD. SILVER SPRING, MD. 


20 

30 


40 

50 

60 


LOGICAL* 1 1128(128,128) 

(M M E N s\oN I S YM (16,3), NSYMI16), IOUT(128,3) 

EQUIVALENCE (IS,LS) , ln/111 . 

DATA NS/16/, NL/3/, IB/1H /, IP/1H+/ 

DATA I S YM/ 

1 1H , 1H. , 1H, , 1H = , 1HU, 1HX , 1H# , 1H8 , 1HM, 1H0, 1H0, 1H0, 1H0, 1H0, 1H0, 

2 1H ,1H 1H IlH IlH IlH , 1H , 1H , 1H , 1H , 1H= , 1H I , 1 H+ , 1 H* , 1 H* , 

3 1H IlH IlH , 1 H , 1 H , 1H , 1H ,1H ,1H ,1H , 1H , 1H ,1H ,1H-,1HX, 
DATA NS YM /l, 1,1, 1,1, 1,1, 1,1, 2, 2, 2, 2, 3, 3, 3/ 

DO 60 J= 1 , 128 
I PMA X= 1 
DO 30 1=1,128 
LS= 1128(1 ,J) 

I F ( IS.LT.l) IS=1 
IF(IS.GT.NS) I S =NS 
DO 20 I L = 1 » N L 
I OUT ( I , IL )= ISYM ( IS, IL ) 

I E ( NSYM ( IS) .GT.IPMAX) IPMAX=NSYM (IS) 

CONTINUE 
I CC= I B 

DO 50 I L= 1 , IPMAX 

WRITE! NGRA, 40) I CC , ( IOUT ( I , I L ) ,1 = 1,128) 

FORMAT ( A 1 , 2 X , 1 2 8 A 1 ) 

ICC=IP 

CONTINUE 

RETURN 

END 


1 HO , 
1H*, 
1HM/ 


E-28 



REFERENCES 


1. Computer Sciences Corporation, 6002-1, Organizational Structures for 
Constant Resolution Earth Data Bases , (prepared for the Environmental 
Prediction Research Facility, Monterey, California), P. R. Beaudet, 

F. K. Chan, and L. Goldshlak, November 1973 

2. Goldstein, H. , Classical Mechanics . Reading, Massachusetts: 

Wesley Publishing Co. , June 1959 


R-l 



